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ABSTRACT 


In the 1990’ s dedicated radio telescopes will be launched into Earth orbit and will be 
integrated in the ground-based VLBI networks. A straightforward extension from present 
ground-based VLBI to space is called space VLBI, which uses radio-antennas in space. In 
the simplest version of the space VLBI technique, one station in orbit observes in 
conjunction with a second station on the ground. However, in practice we have a number 
of networks of ground antennas observing the common celestial radio sources 
simultaneously with a conventional VLBI technique. Moreover, joint observations of two 
or more space VLBI satellites will supposedly be performed in the future. Therefore, in 
our research work a combined use of simultaneous space and ground-based VLBI 
observations are considered from the geodetic and geodynamic point of view. 

This investigation studies the feasibility of space VLBI observables for geodesy and 
geodynamics. A brief review of space VLBI systems from the point of view of potential 
geodetic application is given. A selected notational convention is used to jointly treat the 
VLBI observables of different type of baselines within a combined ground/space VLBI 
network. The basic equations of the space VLBI observables appropriate for covariance 
analysis are derived and included. The corresponding equations for the ground-to-ground 
baseline VLBI observables are also given for a comparison. The simplified expression of 
the mathematical models for both space VLBI observables (time delay and delay rate) 
include the ground station coordinates, the satellite orbital elements, the earth rotation 
parameters, the radio source coordinates, and clock parameters. The observation equations 
with these parameters have been examined in order to determine which of them are 
separable or nonseparable. Singularity problems arising from coordinate system definition 
and critical configuration are studied. Linear dependencies between partials are analytically 

derived. 

The mathematical models for ground-space baseline VLBI observables have been 
tested with simulation data in the frame of some numerical experiments. Singularity due to 
datum defect is confirmed. 

Recommendations are given for future research work. 
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1 . INTRODUCTION 


Two dedicated space VLBI projects are currently in preparation to launch one or 
more VLBI radio telescopes in orbit. The first space VLBI mission that will be 
implemented in the near future is RADIOASTRON in the Soviet Union. It is already an 
approved and funded mission. The satellite will be launched in 1993. Its trajectory is 
planned to be 3,000 km in perigee and 69,000 km in apogee. The satellite will carry a 
10 m antenna. The second project is a Japanese orbiting VLBI mission called VSOP. The 
expected launch date of the satellite is 1995. Its trajectory is planned to be 1,000 km in 
perigee and 20,000 km in apogee. Both projects are now in progress. Most probably, a 
combined use of RADIOASTRON and VSOP observations will be performed in the future. 
There is a Western European mission with NASA participation in Phase A Study at ESA 
called QUASAT. Although this project was not approved by ESA at the selection round of 
October 1988, a second generation space VLBI mission called IVS has been proposed to 
ESA at the end of 1989. 

The feasibility of and potential for using a dedicated VLBI observatory in space has 
already been demonstrated successfully by the NASA TDRSS satellite with a nearly five 
meter antenna. In 1987 during the XIX General Assembly of IUGG in Vancouver, IAG 
Special the Study Group 2.109 — Application of space VLBI in the Field of Astronomy and 
Geodesy — was established to study the usefulness of space VLBI for astrometry and 
geodesy. 

It is expected that space VLBI will be a reality in the current decade. Orbiting radio 
telescopes will be used to make interferometric observations of extragalactic radio sources 
in conjunction with the major ground-based VLBI arrays in Europe, USA, Australia, Japan 
and the USSR. It is planned to determine the orbits of these radio telescopes in space with 
high accuracy by the missions themselves and possibly using additional tracking systems 
(e.g., GPS, PRARE). The main goals of all space VLBI projects are to carry out 
astrophysical investigations. The current missions are devoted to improve imaging quality 
and angular resolution of compact galactic and extragalactic radio sources. At the same 
time, the radioastronomers are going to plan the scientific goals of space VLBI missions to 
be as wide-ranging as possible. 
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The anticipated use of a space VLBI system poses several important questions with 
respect to its usefulness in geodesy and geodynamics. 

The potential applications of the space VLBI system in these areas are related to the 
connection and unification of reference systems and frames inherent in the technique. 
Since the currently available space geodetic techniques are effectively used for these 
applications by station collocation in the frame of the International Earth Rotation Service 
(IERS) [Boucher et al., 1988] and [Mueller, 1988], one should investigate how the space 
VLBI system could complement these techniques. The capabilities of the new space 
techniques to provide valuable information on the systematic differences between the 
frames of the various Conventional Inertial Systems (CIS) and Conventional Terrestrial 
Systems (CTS) is now widely investigated [Kovalevsky et al., 1989] because of the high 
importance of the connection and unification of all reference frames in geodesy and 
geodynamics. In principle, a space VLBI system offers an opportunity to connect the 
reference frames of the CTS and two types of CIS: 1) Radio Source-CIS, and 2) Dynamic 
(that is Satellite Orbit)-CIS inherent in a space VLBI system. The establishment of the 
relationship between them would be of scientific interest. Note that the currently adopted 
Conventional Celestial Reference System (CCRS) is based on radio source positions. 

Another important problem area that needs investigation is that related to the 
usefulness of space VLBI data for gravity field determination. The present gravity field 
determinations using artificial satellites are based mainly on range and range-rate 
observations which contain no direct directional information [Marsh et al., 1988] and 
[Rapp, 1989]. Therefore, they are only indirectly linked to the inertial reference frame. 
One should investigate how the space VLBI observation containing directional information 
could be applied to improve the gravity field determination. The space VLBI missions 
offer and provide new types of satellite observables (VLBI time delay, delay rate, and 
differential VLBI tracking data) with high accuracy for these potential applications as well. 

Earlier investigations on the potential applications of space VLBI for astrometry, 
geodesy and geodynamics include Adam (1989), Adam and Mueller (1989), Bartel (1989), 
Fejes et al., (1986, 1987 and 1989a, b), Kawaguchi (1989), Koyama (1989), and 
Takahashi (1989). 

On the basis of the background mentioned above, it was considered appropriate and 
worthwhile to pursue the present investigation. The aim of which is to briefly review the 
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space VLBI systems in progress emphasizing the geodetic aspects, to derive mathematical 
models for both space VLBI observables (time delay and delay rate) suitable for least 
squares covariance analysis of the parameters of geodetic interest, and to explore 
estimability problems inherent in the space VLBI technique. An additional important aim is 
to numerically check the mathematical models by test computations. 

Consequently, the whole work is divided into five chapters, each of which treats 
aspects of the geodetic applicability of space VLBI. Chapter 2 gives a brief review of space 
VLBI missions in progress. Chapter 3 derives the mathematical models for space VLBI 
observables and studies the estimability of geodetic parameters involved in the models. 
Chapter 4 describes the theoretical basis of the datum problem of a space VLBI network 
within singular Gauss-Markov model and shows the results of test computations. 
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2. BRIEF REVIEW OF SPACE VLBI SYSTEMS 


2.1 Introduction 

The application of the VLBI technique to the investigation of the nuclei of galaxies 
and quasars has allowed the mapping of their angular radio brightness distributions and the 
studying of the kinematics of their cores. Substantial improvements in the resolving power 
of VLBI might enable the proper motion of a galaxy or a quasar to be measured. 
Microarcsecond angular resolution will allow community of radio astronomers to study 
several phenomena, see, e.g., [Bartel et al., 1986 and 1988], [Kardashev and Slysh, 
1988], Annual parallax, proper motion and position determinations (quasar cores, nearby 
galaxy cores, galactic centers, pulsars and radio stars) at microarcsecond (ftas) level are 
also very important for geodesy and geodynamics. 

There are two ways of increasing the angular resolution: (1) use of short radio 
wavelengths (mm-VLBI), and (2) use of radiotelescopes in space (space VLBI). The first 
encounters problems caused by atmospheric phase fluctuations and technical difficulties. 
The second is more straightforward and has no limitations in principle to the baseline 
length, see [Burke, 1983], [Kardashev and Slysh, 1988], [Preston, 1983], [Sagdeev, 
1984] and [Schilizzi et al., 1984], 

A very straightforward extension from present ground-based VLBI is called space 
VLBI, which uses radio antennas in space. Moreover, radioastronomers need multiple 
elements in space (introducing space-space VLBI) in order to go to higher angular 
resolution and/or higher dynamic range mapping. One of such conceptions is suggested by 
Hirabayashi (1989b) for two orbiting VSOP observatory with modestly near earth orbits. 

2.2 Basic Concept of Space VLBI 

The usual geodetic VLBI system shown in Fig. la [Counselman and Shapiro, 1978] 
consists of an array of at least two antennas that observe the same radio source 
simultaneously. A direct electrical connection is not maintained between the antennas, thus 
allowing them to be separated by thousands of kilometers. At each antenna station, the 
radio interferometry signal (RF) received from the observed radio source is converted to a 
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lower, “intermediate” frequency (IF) by mixing with a local oscillator (LO) signal. At each 
site, the IF signal is tape recorded with a reference time base derived from the same 
frequency standard. Tapes recorded simultaneously at the antenna sites are later collected 


ANTENNA 1 ANTENNA 2 



Fig. la Ground-based geodetic VLBI system 

and played back together at a central processing station where the reproduced signals are 
cross-correlated to obtain the basic VLBI observables. 

The space VLBI system does not differ conceptually from the case of a ground-based 
array of radio telescopes. A schematic diagram of the space/ground VLBI system is 
represented in Fig. lb [Schilizzi et al., 1984] and [Schilizzi, 1988a]. The space-borne 
antenna will observe the same radio sources in conjunction with networks of antennas, and 
relay the received signals via a digital or analogue link directly to telemetry stations on the 
ground. A phase/frequency reference for the antenna in space will be based on hydrogen 
maser oscillators on the ground and relayed directly to the satellite from the telemetry 
stations in turn (phase transfer). The stability required for this phase transfer is very high 
(about lxlO' 14 ). After transmission to the ground, the IF data will be recorded on VLBI 
magnetic tapes on the ground in exactly the same way as for the ground based elements of 
the array. The tapes will then be brought together with tapes from the ground VLBI arrays 
at a central processing station for cross-correlation and image processing. After cross- 
correlation and calibration, the obtained basic VLBI data will be used for scientific 
investigations. 
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Fig. lb A schematic diagram of a Space/ground VLBI system 

All communication with the space VLBI antenna will be through one or more 
telemetry/control (T/C) stations in the network. The two way or multiple way coherent 
link(s) by telemetry stations provide range, range-rate, and phase data which can be used 
for orbit determination. 

2.3 TDRSS-OVLBI Demonstration Experiment 

The feasibility of and potential for using a dedicated VLBI observatory in space has 
already been demonstrated (see [Levy et al., 1986, 1987 and 1989], [Linfield et al., 1988 
and 1989], [Hirabayashi, 1988], and [Nishimura and Hirabayashi, 1988]. A 
communications satellite in the Tracking and Data Relay Satellite System (TDRSS) of 
NASA was successfully used for the first space VLBI test experiment in 1986-87. The 
idea of using existing TDRSS for demonstrating space VLBI has been proposed by 
scientists at the Jet Propulsion Laboratory (JPL). The TDRS has two identical high-gain 
antennas with a 4.9 m diameter for S- and K-bands (see Fig. 2). The TDRSE which is 
situated at 41° W in the geostationary orbit was used in combination with ground 64 m 
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diameter antennas in Tidbinbilla, Australia, and in Usuda, Japan. It was the most suitable 
existing satellite for orbiting VLBI (OVLBI) because of the design of its local oscillator 
chain and its sensitivity. The test experiment called TDRSS-OVLBI experiment was done 
with the observing frequency of 13 cm (2.3 GHz). The received signal from radio sources 
at TDRSS is down-linked to the ground, frequency converted, and recorded on magnetic 
tapes. The reference signal in turn is uplinked to the satellite to phase lock the receivers on 
board to maintain coherence of the interferometer. Fig. 3 shows the concept of the 
TDRSS-OVLBI demonstration experiment. 

The TDRSS-OVLBI experiment was performed in two phases (July-August, 1986 
and January, 1987). For the observation, radio sources in the Southern Hemisphere 
ranged in declination from 0°to -30° were selected due to the limits of the TDRSE satellite. 
In the July-August 1986 session, fringes were successfully obtained among three stations 
(TDRSE, Usuda and Tidbinbilla) for three quasars. The maximum projected baseline 
obtained was 17,800 km (about 1.4 times Earth’s diameter) by the TDRSS-Usuda 
baseline. In January, 1987, additional observations were made using the same ground 
antennas. Of the 25 observed sources, 23 were detected on orbiter-ground baselines, with 
baseline lengths as large as 2.16 Earth diameters. Also, the radio source structures were 
estimated from the January, 1987, session. Note that in 1988 an experiment at 2.3 GHz 
(13 cm) and 15 GHz (2 cm) frequencies was also successful. 

Tracking information for TDRSE is obtained through triangulation with ground 
transponders located at White Sands Ground Terminal (WSGT) in New Mexico (32° N 
latitude, 107°W longitude) and Ascension Island (8° S latitude, 15° W longitude) (see 
[Levy et al., 1989], The forward link uses a pseudo-noise code. The return link uses a 
different but synchronized pseudo-noise code modulation. By comparing epochs of the 
pseudo-noise code (all measurements are referenced to the WSGT cesium standard), the 
range is obtained. The Doppler frequency is also extracted and counted to obtain the range 
rate. The ranging and Doppler tracking data were processed by K.B. Blaney at Goddard 
Space Flight Center (GSFC) for orbit determination by a differential correction technique 
(ibid.). The calculation procedure is a least-squares fitting of weighted range and Doppler 
residuals to a corrected orbit. A data span (and solution interval) of 30-34 hr was used. 
The solution yielded three position components, their associated velocity components, and 
their time derivatives at 60 s intervals. 
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As a summary, the test of the technical concept and the probe of the scientific 
potential for space VLBI were successful. The TDRSS-OVLBI experiment shows that an 
accurate knowledge of the spacecraft orbit is essential in order to achieve the coherence of 
the measured data with high sensitivity and calibration accuracy. 

2.4 Future Space VLBI Missions 

Projects to place one or more telescopes in space for VLBI purposes have been 
proposed or planned by several space agencies: 

• QUASAT (Quasar Satellite) by ESA and NASA, see, e.g., [Schilizzi et al., 1984] 
and [Schilizzi, 1988a] 

• RADIOASTRON (Radio Astronomical Satellite) by The Space Research Institute 
(IKI) of the Academy of Sciences in the USSR, see, e.g., [Sagdeev, 1984], 
[Kardashev and Slysh, 1988] and [Kardashev, 1989] 

• VSOP (VLBI Space Observatory Program) by The Institute for Space and 
Astronautical Science (ISAS) in Japan, see, e.g., [Hirabayashi, 1984 and 1988] and 
[Nishimura and Hirabayashi, 1988], 

• I VS (International VLBI Satellite), a proposal to ESA to follow RADIOASTRON and 
VSOP [Schilizzi et al., 1989]. 

The orbital parameters of the planned space VLBI satellites are given in Table 1. 

2.4.1 The RADIOASTRON Mission. The first space VLBI mission that will 
be implemented in the near future is RADIOASTRON. It is an approved and funded 
mission in the Soviet Union. The design study is in progress with significant contributions 
from abroad. The RADIOASTRON satellite will carry a 10 m antenna, and its trajectory is 
planned to be 3,000 km in perigee and 69,000 km in apogee. The launching is expected 
around the end of 1993. 

Fig. 3 is an overall schematic diagram of the Radioastron system. The satellite radio 
telescope (Fig. 4) will be a prime focus system with a deployable parabolic reflector 10 m 
in diameter. The reflector has a fixed inner part of 3 m diameter and 24 unfoldable panels. 
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Table 1. Comparison of the Orbital Parameters of QUASAT, 
RADIO ASTRON and VSOP Satellites. 



QUASAT 

RADIOASTRON 

VSOP 

h a (km) 

36,000 & 22,000 

69,000 

20,000 

hp (km) 

5,000 & 5,000 

3,000 

1,000 

T (hrs) 

12.2 & 7.75 

24 

6.06 

i(°) 

30 

65 fixed 

46.4 or 31.0 

o>(°) 

variable 

285-315 
to be chosen 

variable 

e 


0.82 

0.379 

Q. 


to be chosen 

variable 


h a = apogee height, h p = perigee height, T = period, i = inclination 

to = argument of the periapsis, e = eccentricity, O = longitude of the ascending node 


The reflector is made of reinforced carbon fiber and has an rms surface accuracy of 
0.5 mm. Four dual circular polarization receivers will be provided on the satellite. The 
spacecraft operation time will be from two to five years. 

The down-link to the ground station is at X-band (8.2 GHz) with a two-channel 
transmitter and high-gain parabolic antenna. The ground-based station of the 
RADIOASTRON mission in Suffa, USSR consists of one transmitting (command and S- 
band local oscillator up-link) antenna and two receiving antennas (S-band local oscillator 
down-link and X-band high-data rate down-link). The 70 m antenna will be used as a 
radio telescope to receive signals from a radio source. Both X-band data and 70 m radio 
telescope outputs will be recorded on two tape recorders in the VLBA format. 
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2.4.2 The VSOP Mission. The Japanese VSOP satellite will carry a 10-m 
antenna, and its trajectory is 1,000 km in perigee and 20,000 km in apogee. This orbit 
provides good imaging capability of VSOP with the combination of the telescope arrays on 
the Earth. Fig. 6 is a rough sketch of the VSOP satellite. The surface accuracy of the 10m 
antenna will be better than 0.5mm RMS to be operable up to 22 GHz (wavelength of 1.3 
cm) reception. The pointing accuracy of the high gain antenna is planned to be 0.01° and 
thus a high precision star sensor will be employed. The local oscillators on-board the 
VSOP satellite must be phase-locked to the frequency standard on the ground. VSOP will 
have cooled amplifiers, giving more than one order of magnitude better sensitivity than 

TDRS. 

The position and the velocity of the satellite will be tracked much more precisely than 
other satellites by a sophisticated tracking scheme. VSOP is now in progress and the 
expected launch date is 1995. The space VLBI program VSOP itself has been informally 
approved by the Japanese government very recently [Hirabayashi, personal 
communication, May 30, 1989]. More details about VSOP mission can be found in 
[Hirabayashi, 1984; 1988; 1989 a, b, c] and [Nishimura and Hirabayashi, 1988]. 



Fig. 6 The VSOP satellite [Nishimura and Hirabayashi, 1988]. 



2.4.3 The QUASAT Mission. The proposed QUASAT mission is a free- 
flying spacecraft carrying a 15 m radio antenna (Fig. 7). Two orbital situations have been 
considered for the mission. One of the orbital situations is featured by the first operational 
elliptic orbit with perigee of about 5,000 km and apogee of about 36,000 km. It was 
planned that the spacecraft’s apogee be lowered later to 22,000 km in order to fulfill the 
various scientific requirements. The orbiting radio telescope was planned to make 
interferometric observations of radio sources in conjunction with the major ground-based 
VLBI arrays in Europe, USA, Australia and the USSR. The mission design lifetime was 
planned to be two years, but an operational lifetime of five years was expected. More 
details about QUASAT mission can be found in, e.g., [Proc. of Workshop on QUASAT: 
1984], [Schilizzi et al., 1984], [Frisk et al., 1988] and [Schilizzi, 1988a], 



Fig. 7 The QUASAT satellite [Schilizzi, 1988a]. 


2.4.4 The IVS Mission. The IVS mission is conceived to be a major radio 
telescope in space funded by the principal space agencies. Orbiting the Earth, the 25 m 
diameter telescope will provide high quality images of radio sources to wavelengths 
spanning the radio band from meters to millimeters will resolutions as high as 25 micro 
arcseconds. It is proposed that IVS be placed initially in an orbit whose apogee height is in 
the range 20,000 km (VSOP-like) to 40,000 km (QUASAT-like), but whose inclination is 
higher than either VSOP or QUASAT. Possible orbital parameters are therefore: apogee 
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height 20,000-40,000 km, perigee height about 5,000 km and inclination about 60 degrees. 
Maximum baseline lengths of 30,000-50,000 km can be generated in this manner, three to 
five times longer than the maximum available on Earth. IYS is a second generation mission 
well matching the ground VLBI arrays in both sensitivity and wavelength coverage. The 
goals of the IVS mission will require the position of the satellite to be known to 30 cm or 
better at all times [Schilizzi et al., 1989]. 

2.4.5 Joint Space VLBI Observations. The space VLBI missions are 
expected to be operational in the next decade. However, QUASAT was not approved by 
ESA at the end of 1988. Therefore, the Soviet RADIOASTRON and the Japanese VSOP 
will be the first two dedicated space VLBI missions in the 1990’s. Supposedly, QUASAT 
will again be considered by ESA to realize. These missions will have many improvements 
compared to the TDRSS which was used for the first space VLBI test experiment. The u-v 
coverage in the TDRSS-OVLBI experiment was not satisfactory because of the 
geostationary nature of the satellite and of the small number of ground radio telescopes 
which joined the experiment. 

The RADIOASTRON and VSOP orbits are complementary in the u-v coverage they 
will provide. A simulation study by Murphy (1989) shows that joint space-VLBI 
observations can produce better results in imaging processing. Therefore, further 
simulation tests will be carried out to determine the optimal orbital elements of 
RADIOASTRON for possible simultaneous observation with VSOP. Their combined use 
is most probable and desirable. It is expected that a combined use of 
RADIOASTRON and VSOP observations will be performed in order to have 
much improved images. 

Radioastronomers need multiple elements in space for higher angular resolution 
and/or higher dynamic range mapping. There are plans to have multiple elements in space 
within both missions introducing space-space VLBI. Hirabayashi (1989b) presented a plan 
for two orbiting VSOP observatory with modestly near earth orbits. VSOP mission with 
two space elements might be very appropriate particularly for geodetic and geodynamical 
purposes as well. 

The coordination for international collaboration is a critical item for both projects — 
RADIOASTRON and VSOP. The design study of both missions is in progress with 
significant contributions from abroad. For space VLBI to become operational, 
coordination of observing, management of ground telescopes and tracking operations, data 
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transfer and processing is needed. Therefore, international workshops are organized. The 
first was to discuss QUASAT in June 18-22, 1984 in Austria fProc. of Workshop in 
QUASAT 1984], The International Radioastron Review Meetings have been held twice a 
year. The eighth meeting was held in Green Bank, West Virginia, USA [Memorandum of 
the 8th RADIOASTRON Meeting 1989]. The ninth one was in Tashkent, USSR in 
November, 1989. Domestic VSOP meetings have been held in Japan [Proc. of a Domestic 
Workshop on “The Research of Space VLBI” 19891. An international symposium was 
organized also in Japan December 5-7, 1989 in order to discuss worldwide the VSOP 
mission. An other international meeting “Space VLBI — The Missions and the Science 
organized by Committee on Space Research (COSPAR) in The Hague, Netherlands, July 
2-3, 1990 reviewed all aspects of space VLBI, with special emphasis on the two approved 
missions RADIOASTRON and VSOP. 

2.5 The Ground-Based VLBI Networks 

The ground-based arrays are the other major element of the space VLBI missions. 
All missions assume an orbiting observatory to be coherently connected to the 
radiotelescopes distributed widely on the earth. They require a global network of ground 
antennas to be able to form subnets. Therefore, the involvement of all possible ground 
antennas is important not only for radioastronomical purposes but for application of them in 
the field of geodesy and geodynamics as well. 

There are two main networks in operation — the European VLBI Network (EVN) 
and the U.S. Very Long Baseline Array (VLB A), both of which provide observing 
opportunities on regional as well as international scales, see [Romney, 1988) and 
[Schilizzi, 1988b]. Networks of both types, tracking/link and observing, respectively in 
the USSR and in Japan, are very important from the viewpoint of their space VLBI 
missions. They will constitute the main support to the orbiting elements. However, 
other networks in the Southern Hemisphere (e.g., the Australia Telescope (AT) and 
Australian VLBI antennas), in Asia (e.g., China, India) and in Canada are also likely to 
participate. Stations of these networks will be the co-observing stations for the space VLBI 
missions. 

Note that VLB A construction is now in progress with expected completion of 1992. 
EVN and AT constructions are also in progress. A construction of dedicated VLBI 
network “QUASAR” of six radiotelescopes in the Soviet Union with possible extensions 
by other telescopes from abroad (China, India) has been started very recently | Finkelstein 
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and Yatzkiv, 1989]. These VLBI networks (and correlators) are expected to work jointly 
with space VLBI missions. 

Tracking networks like the NASA/JPL Deep Space Network (DSN) will collaborate 
on precise tracking of both RADIOASTRON and VSOP satellites. The different networks 
of VLBI antennas are shown in Figs. 8 and 9, which reflect the strong international 
character of this research. The observing schedule, data processing, the managements of 
ground telescopes and tracking operations need to be discussed and organized worldwide. 

2.6 Orbital Tracking of Space VLBI Satellites 

The astrometric, geodetic and geodynamic potential of the space VLBI missions is 
strongly dependent on orbit determination accuracy of the satellites. Therefore, high 
accurate knowledge of the satellite orbit is very important for the space VLBI missions. The 
following tracking techniques are recommended for precise orbit determination of space 
VLBI satellites, see [Fejes, 1989] memorandum of the 8th RADIOASTRON Review 
Meeting 1989 [Tang, 1984], 

The first category which requires no systems modifications or additional equipment: 

1 . Two-way or multiple-way coherent link(s) by telemetry stations (range, range- 
rate, phase). A good telemetry coverage is essential. 

2 . VLBI delay and delay-rate as tracking data. 

3 . Differential VLBI (A VLBI) using extended (worldwide) network of ground radio 
telescopes during selected observing programs. 

Results of a computational simulation study by Konopliv (1989) show that it might 
be possible to reach meter-level position accuracy in orbit determination of 
RADIOASTRON satellite with range and AVLBI measurements by NASA/JPL DSN 
stations. Borza et al., (1990) show by simulations that the dominant part of the orbital 
errors of space VLBI satellites originates from errors in the solar reflectivity models. All 
other effects can be modeled so accurately that they only yield small or negligible errors. 
According to the test computations carried out by Fejes et al., (1989b), the space VLBI 
time delay and delay-rate observables can contribute significantly to orbit accuracy 
improvement for the missions. It is in a good agreement with the statement of Tang 
(1984). 
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The second category which requires system modification and additional ground and 
space equipment: 

1 . Microwave techniques: (a) One-way satellite-based systems (GPS, GLONASS); 
and (b) two-way satellite-based systems (PRARE) 

2. Laser ranging 

3 . Micro-accelerometer on board 

Note that the inclusion of PRARE System for Radioastron is considered. It is 
planned that the VSOP satellite will carry a Global Positioning System (GPS) receiver on 
board to provide tracking data for orbit determination. A set of comer-cube reflectors on 
board of a space VLBI satellite (e.g., VSOP) would increase the potential of space VLBI in 
geodesy and geodynamics. 
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3 . MATHEMATICAL MODELS FOR SPACE VLBI OBSERVABLES 


3.1 Introduction 

Since there are no conceptual differences between ground-based and space-based 
VLBI, existing formalism of the ground-based geodetic VLBI carries over easily to the 
space-based case. In ground-based VLBI the necessary relative motion of the stations is 
provided by the Earth’s rotation. In space-based VLBI the relative motion of the antennas 
is provided by both Earth’s rotation and orbital motion of the space VLBI satellite around 
the Earth. 

In the simplest version of space VLBI, one station in orbit observes in conjunction 
with a second station on the ground. However, in practice we have a number of networks 
of ground antennas observing the common celestial radio sources simultaneously with a 
conventional VLBI technique. Moreover, joint observations of two or perhaps more space 
VLBI satellites will supposedly be performed in the future. Therefore, a combined use of 
simultaneous space and ground VLBI observations should be considered. 

3.2 Basic Observables and Their Models 

We have selected notational conventions which allow us to treat the VLBI 
observables of the different types of baselines within a combined ground/space VLBI 
network uniquely. A space VLBI network is formed by ground VLBI stations Pj,Pj and 
orbiting VLBI satellites’ positions S 1 and S J . The lower index with small letters (i,j) refers 
to ground points while the upper index with capital letters (I,J) to satellite points at high 
altitude. Therefore, in our conventions we can distinguish a ground-to-ground baseline 
between two ground VLBI stations Pj and Pj by Bjj , a ground-to-space baseline between 
a ground VLBI station Pj and an orbiting VLBI telescope S J by B{ , and a space-to-space 
baseline between two orbiting VLBI telescopes S 1 and S J by B u . We use these notations 
to distinguish the corresponding VLBI observables, i.e., time delay and delay late as well. 

The basic geometry for a typical baseline in the ground-based VLBI and space-based 
VLBI or simply space VLBI is shown in Figs. 10-12. Generally, a certain segment of a 
wavefront from the natural radio source will arrive at one site before it arrives at another 
(Figs. 10-12). This time delay is the basic observable of VLBI. Delay observations are 
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composed of three components. The first component is due to the geometry of the station 
locations and source location, as illustrated in Figs. 10- 12. This component contains all of 
the information about the geodetic and astrometric parameters that are of interest. The 
second component is due to instrumental effects, principally clock errors. The third 
component is due to effects of the propagation medium (e.g., the atmosphere and the 
ionosphere). This latter component is not considered here. 

In the ground-based VLBI network, two ground radio telescopes Pj and Pj form a 
ground-to-ground baseline By (Fig. 10). The signal from the ^ th radio source arrives at the 
k* epoch of observation to the ground VLBI stations with time delay tyk/ which can be 
written as 


'Cijw = ~ djjw = Bjj(tk) cos y, , 
or in vector notation, the dot denoting the inner product, 

Ljk/ 1 = — ~ By (tk) • tl 


where e* is the unit vector in the direction of the ^ th radio source, By is the vector 
separation of the receiving stations Pj and Pj, and c is the speed of light. The minus sign is 
introduced to follow convention. The changing geometrical configuration gives rise to the 
delay rates TyM 


^ijk£ 


1 

C 


dBy (t k ) 

dT~' e * 


(3) 


assuming that e* = 0 . 


The third VLBI observable is phase delay which is related to time delay (2) 

by 

where to is the angular frequency (2nf) of the received radio signal. The phase delay is the 
difference in phase of the signals received at each site. 

Note that time (group) delay Tyw is the derivative of (fringe) phase 0ijk/ with 
respect to angular frequency. The time delay contains the same information as phase for 
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astrometric and geodetic purposes but is inherendy less accurate than phase. On the other 
hand, it can be estimated unambiguously. The delay observable yields a full baseline 
solution and, therefore, plays the most important role in geodetic VLBI. The delay rate is 
not often used in ground-based VLBI. However, the role of the rate observable will be 
more important in space VLBI. 

In the following the speed of light will be set to unity so that the time delay and 
delay rate will be expressed in units of distance and its time variation, respectively. 
Assume that they have been corrected for discrepancies between modeled and true distance 
(or distance variation), e.g., atmospheric errors, relativistic effects and abberation. The 

time delay is now 

dijw = c Tyu = - By(tk) • e/ , (5) 


and the time delay rate is 


. . dB;j(tk) 

dijk* = c ^ijk t ~ ^ 


• Zl ■ 


( 6 ) 


Adding a two term polynomial, whose coefficients ACoy and ACi- correspond to a relative 
offset and rate, respectively, between the two clocks at the ends of the baseline, the 
geometric time delay (expressed in units of length) can be modelled as 

djjw = - By(tk) -e< + c [AC 0ij + AC 1;j (t k - t G )] , (7) 


where to is the initial epoch of observation. The time delay rate is then 


dBi,(tk) A/ ~, 

d ijk / = + cACiy • 


( 8 ) 


The clock parameters AC 0ij , AC ly are nuisance parameters, included here to make the 
model more realistic. 


In space VLB! network, an Earth-orbiting VLBI telescope S 1 with two ground 
VLBI stations Pi and Pj at a moment t k forms two ground-to-space baselines B }(t k ) , Bj(tk) 
in addition to the ground-to-ground baseline Bij(tk) (Fig. 11). The corresponding time 
delays At , tfk/ and delay rates At , tfkt , including the clock parameters, can be 


22 



expressed by similar observation equations which are for the ground-to- space baseline 
Bj(tk) as follows: 


4 rc^r- Bj(tk> c/ + c [AC l o r} + AC} rj (t k - to)] 


and 


d]k/ = c Xjw = 



■ei + c ACj rj 


(9) 

( 10 ) 


where the delay rate (tfu) plays a more important role than at the ground-based VLBI 
network due to the larger relative motions between a ground and an orbiting antenna. 

The lower index r of the clock parameters in Eqs. (9) and (10) refers to the 
reference clock at the telemetry/control (T/C) ground station. There will be no clock on- 
board of the orbiting radio telescope (see Chapter 2.2). 

Let us suppose that there will be two orbiting space VLBI satellites in operation in 
the future (e.g., RADIOASTRON and VSOP, IVS may be later in addition, or two VSOP 
satellites as a system in themselves) having simultaneous observations for the same 
extragalactic radio source. Two orbiting VLBI satellites, S 1 and S*, at a moment t k form a 
space-to-space baseline B IJ (t k ) (Fig. 12) . Then, including the clock parameters, the 
corresponding time delay and delay rate are modelled 

<4* = c xfi = - B IJ (tk) • e, + c[ACo rq + AC} J rq (t k - to)] (U) 


dw = c iu = 


dB IJ (t k ) 

dt 


t t + cAC? rq , 


( 12 ) 


where the lower indexes r and q refer to the clocks at telemetry/control ground stations r 
and g, respectively for the satellites S 1 and S J . 

By examining the geometric part of the time delay observable models we recognize 
that the time delays will change if the length of the baseline changes, or if the orientation of 
the baseline changes with respect to the plane of the wavefront, or if both change 
simultaneously. In the case of a ground-to-ground baseline, the time delay changes due to 
the variation of baseline orientation with respect to the plane of the wavefront. In the 
ground-based VLBI, the length of a baseline is supposed to be unchanged for a short 
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period (the geodynamical effects, e.g., tides, can be accounted for). It is evident that time 
delay (and delay rate) observations can be used to recover baseline length and baseline 
orientation with respect to the radio source direction or recover the source orientation with 
respect to the baseline direction. The first two applications are of geodetic interest while the 
last one is of astrometric. 

In the space VLBI the time delay changes because of the continuous changes of 
baseline length and of baseline orientation with respect to the plane of the wavefront. Both 
changes are due to the continuous relative motion between both ends of the baseline. For 
the case of a ground-to-space baseline, the relative motion between both ends of the 
baseline is a result of the Earth’s rotation and satellite motion. Changes in both length and 
orientation of a space-to space baseline, result from the relative motion of satellites with 
respect to each other. Since the baselines in space VLBI do not repeat themselves, their 
length and orientation therefore can not be recovered from the usual VLBI technique itself. 
Some of the parameters (e.g., satellite’s coordinates and/or source coordinates) should be 
determined by some other means. 


3.3 Datum Content of Space VLBI Observables 


Following an evaluation by Wells et al., (1987) for GPS datum definition, in this 
chapter we describe datum content of the space VLBI time delay observables. For this 
purpose, consider a space VLBI network of points divided into satellite points with an 
orbiting radio telescope S^ and a ground-based VLBI station Pj. A single time delay 
observation (ignoring clock parameters, and expressing the time delay in units of distance) 
is given by 


djk/ — ' 


TT 

RkXj X{ , 


in matrix form, where Xj contains the Earth-fixed coordinates of the station Pj, X k contains 
the coordinates of satellite S 1 in a true-of-date geocentric inertial coordinate frame at epoch 
t k , R k is a transformation matrix (transpose of the matrix S), see Eqs. (7-11) in [Moritz 
and Mueller, 1987; p. 417], the dot denotes the scalar product and the T is the transpose 
sign. 


Translating the datum by AX we have 

R^Xp R k Xj + AX 
x{ = Xjj + AX 
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(13a) 


_ r _ it 

djW = — _R k Xj — X k 

= -.(R k Xj + AX)-(Xj[ + AX)] 

= - R k Xj - x]y ti 

= d]k<f 

that is, time delays are invariant with respect to translation. 

Rotating the datum, we have 

RkXj = Ro(RkXj) 

X k = RoX^ 
ti = R 0 t e 

where R<, represents a proper orthogonal transformation (R„ R 0 = I) . Then 

~i r. — ~tTL 

djk/ = -[R k Xj-XkJ ti 

1 T 

= — R 0 (R k Xj) — R 0 Xj Ro^/ 

= - R k Xj-Xt] RjR 0 e, 

= -RjR 0 [R k Xj-xf e, 

= - [R k Xj - xjy t e 

=d h 

that is, time delays are invariant with respect to rotation. 

Changing the datum scale by s , we have 

RkXj = s R k Xj 

Xk = sXj 
e* =se ; 

Then 
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( 15 ) 


djk/ = - [RkXj - Xkj 1 ^* 

= - [sR k Xj - sXjj set 
= — s[R k Xj — xf s®/ 

= — s^RkXj — xtf ^ 

= s 2 dj k/ 

that is, time delays are not invariant with respect to scale change. 

For the sake of completeness and a comparison, the corresponding equations for 
the ground-based VLBI time delay observable can be derived in order to show its 
translational and rotational invariance. In the following only the translational invariance 
will be deduced pointing at some difference in the sequel. Consider a network of ground- 
based VLBI stations Pi and Pj. A single time delay observation (ignoring again clock 
parameters and expressing the time delay in units of distance) is given here in matrix form 

dijk/ = -(Xj-Xi) T R k e,, 

where Xj and Xj is the ith and jth station within a ground-based geodetic VLBI network, 
R k is a transformation matrix, see Eq. (19). 

Translating the datum by AX we have 

Xj = Xj + AX 
Xj = Xi + AX . 


Then 

dijk/ = - (Xj - Xi) T R k e* 

= - [(Xj + AX) - (Xi + AX)f R k e< 

= - (Xj - Xi) T R k 

= diju (13b) 

that is, ground-based VLBI time delays are invariant with respect to translation. A 
comparison of Eqs. (13a) and (13b) is discussed in section 3.5.1. 
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Therefore, time delay observations cannot provide datum origin and orientation. 
Consequently, the basic formulations are independent of the reference frame used. Since 
the geometric observables are the dot product of two vectors (B(t) and e/, dB(t)/dt, and e;, 
respectively) it is evident that they do not depend on the location and orientation of the 
coordinate frame in which the components of these vectors refer. 

Eqs. (7-12) can be applied for covariance analysis either in terrestrial frame or in the 
true celestial frame. In the former case, e/ is transformed from the true celestial to the 

terrestrial frame, and in the latter case the baseline vector (By , B* or B IJ ) is transformed 
from the terrestrial to the true celestial frame. On the other hand, the observables are 
affected by the scale of the coordinate frame which is implied by the adopted speed of light. 
Therefore, from observations of both time delay and delay rate it is impossible to recover 
the origin and absolute orientation of the coordinate frame. 

3.4 Mathematical Model for Ground-to-Ground VLBI Observations 

The basic theory of ground-based geodetic VLBI is very well documented by, e.g., 
Whitney (1974), Robertson (1975), Counselmann (1976), Counselmann and Shapiro 
(1978), Ma (1978), Shapiro (1978), Bock (1980), Brouwer (1985), Harvey (1985), 
Campbell (1987), Moritz and Mueller (1987) and Sovers and Fanselow (1987). Therefore, 
only those aspects of the theory which are essential to the research in this report, are briefly 
summarized. 

Since a combined use of ground-based and space VLBI observables will be used in 
the future, a detailed comparison of the corresponding mathematical models is desirable. 
Only the essentially necessary equations for ground-to-ground baseline observables are 
given here in order to make a comparison between the mathematical models for different 
type of baseline observations. 

Therefore, some of the important mathematical equations will be derived in rigorous 
manner. 

3.4.1 Time Delay Model. In the ground-based geodetic VLBI, the mathematical 
model for the time delay observable can basically be written as the inner product of the 
baseline vector rotated from an Earth-fixed system into the true-of-date inertial system and 
the quasar unit vector in the true-of-date inertial system plus polynomial terms of clock 
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parameters. For baseline PiPj observing source t at epoch k the path difference (time delay 
times the speed of light) which the incoming signal must travel after its reception at station 
Pi till its arrival at station Pj (see Fig. 10) can be expressed in a simplified manner suitable 
for least squares covariance analysis as follows (the T is the transpose sign), 

AXjj 1 T [ cos8* cosa* 

d ijki =- AYij R2B)RiH)R 3 (e k ) cos5* sina^ + (16) 

AZij J L sin5 * -I 

+c[AC 0ij + ACi ij (t k -t o )j, 

where 

are the coordinate differences of the baseline PiPj in an Earth-fixed 
system; 

are the true right ascension and declination of the ^th quasar, 
are the components of polar motion that relate the true celestial 
pole to the terrestrial pole; 

is the Greenwich Apparent Sidereal Time (GAST) at epoch t k , that 
is 

(TAI - UTC) - (UTC - UTl)] k + Eq.E ( 17 ) 

6 0 GAST at initial epoch to 

TAI International Atomic Time 

UTC Coordinated Universal Time 

UT1 observed Universal Time corrected for polar motion 
Wd conversion factor from Universal to Sidereal Time 

Eq.E Equation of Equinox 

ACojj, ACijj are the clock offset and drift between the ground stations Pj 
and Pj, 

Ri(<|)) is the rotation matrix for a right handed rotation <(> about the axis i, see 
[Moritz and Mueller, 1987; p. 417]: 


AXij, AYij, AZij 
a l, 5 1 

6k 

e k =e 0 +w d uTi k 
= e 0 + w d [tai - 
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R2(-4) = 


RiC-n) = 


R 3 (0k) = 


cos ^ 

0 

sin E, 

0 

1 

0 

-sin E, 

0 

cos E > 

1 

0 

0 

0 

COS Tj 

-sin ij 

0 

sin Tj 

COS Tj 

cos 0k 

sin 0k 

0 

-sin 0k 

cos 0k 

0 

0 

0 

1 


that is 


R = R 2 U)Ri(-ti) R 3 (0k) = 


(18a) 


(18b) 


(18c) 


cos E, cos 0k-sin E, sin Tj sin 0k cos E s sin 0^+sin ^ sin Tj cos 0k sin E, cos Tj 
-cos Tj sin 0k : cos Tj cos 0k : -sin Tj 

-sin E, cos 0k— cos ^ sin Tj sin 0k -sin E, sin 0k+cos £ sin Tj cos 0k cos £, cos Tj 


Expression of Eq. (16) may be rewritten as 

djjk/ = - AXjj [cos £, cos 8* cos (0k - a e ) - sin E, sin T| cos 5/ cos (0k - oc^) 
+ sin E, cos T) cos 8^] + 

+ AYij [cos T| cos 8/ sin (0k - a^) + sin Tj sin 8/] + 

+ AZjj [sin E, cos 8; cos (0k - a^) + cos E, sin rj cos 8/ sin (0k - a^) 
- cos E, cos Tj sin 8^] + 

+ c ^ACojj + ACijj (tk - to)l 


(19) 


( 20 ) 


Differentiation of Eq. (20) yield the usual form of the adjustment model (error equations), 


d(djjki) — ^ Ap dPp, 

p 

where the index p stands for the unknowns 

[AXij, AYij, AZij t at, 8 t , Tj, k, AC^, ACi-} 


( 21 ) 


( 22 ) 


The parameter k is used to model the departure of the earth’s sidereal rotation from 
uniformity, thus, e.g.. 
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K = UT1 - UTC. ( 23 ) 

The partial derivatives A p are as follows: 

A AXij = - [cos % cos 8* cos (0 k - a ( ) - sin £, sin Tj cos 8; cos (0 k - a*) 

+ sin b, cos Tj cos 8;] (24) 

A A Yjj = cos 'H cos sin (0 k - a*) + sin T) sin 8* (25) 

A A Zjj = sin % cos cos (0 k - ct/) + cos b, sin Tj cos 8/ sin (0 k - <Xi) 

- cos b, cos Tj sin 8^ (26) 

A^ = AXij [sin b, cos 8; cos (0 k - a^)+ cos £, sin ri cos 8/ sin (0 k - a*) 

- cos £, cos T| sin 8/] + 

+ AZy [cos b, cos 8/ cos (0 k - ct^) — sin £, sin Tj cos 8/ sin (0 k — a^) 

+ sin E, cos TJ sin 8*] (27) 

A^ = AXjj [sin % cos T) cos 8^ sin (0 k - a^) + sin b, sin r| sin 8/] - 
- AYy [sin Tj cos 8; sin (0 k - a*) - cos T| sin 8/] + 

+ AZy [cos b, cos ri cos 8; sin (0 k - a^) + cos b, sin Tj sin 8/] (28) 

A k = Wd {AXjj [cos £, cos 8; sin (0 k - a e ) + sin b, sin ri cos 8/: cos (0 k - a*)] + 

+ AYjj cos r| cos 8; cos (0 k - a^) - 

- AZjj [sin b, cos 8* sin (0 k - a*) - cos S, sin Tj cos 8* cos (0 k - a^)] } (29) 

A a( = - {AXjj [cos E, cos hi sin (0 k - a e ) + sin E, sin Tj cos 8; cos (0 k - a^)] + 

+ AYij cos Tj cos 8; cos (0 k - a^) - 

- AZij [sin b, cos 8* sin (0 k - a^) - cos E, sin ri cos 8; cos (0 k - a^)] } (30) 

Ag, = sin hi {AXjj [cos E, cos (0 k - a e ) - sin E, sin r[ cos (0 k - a,)] - 

- AYij cos Tj sin (0 k - a i) - 

- AZij [sin E, cos (0 k - a t ) + cos b, sin Tj sin (0 k - a^)] ) - 

- cos hi [AXij sin ^ cos Tj - AYjj sin Tj + AZij cos b, cos Tj] ( 31 ) 

AACo ij = c (32) 

AACiy = c (tk - to) (33) 
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The set of unknown parameters in Eq. (22) can usually be determined in a standard 
least squares adjustment from the information present in the observables. It is not the 
purpose of this report to describe the adjustment procedure in detail. However, the 
estimability aspect of the adjustment problem will be emphasized. In the present context 
the term simply refers to parameters which can be estimated through the adjustment 
process. The detection of what is estimable and what is not is through the design matrix A 
of the partial derivatives. Therefore, it has particular significance in the following. If there 
exist linear relationships between the columns of the design matrix A, its column rank will 
be deficient (will not be full) and the normal matrix A T PA (P is the weight matrix of the 
observables) will consequently be singular implying that not all of the parameters are 
estimable and the establishment of a new set of the parameters is required. 

The following linear relationships can be detected among the partial derivatives in 
Eqs. (24) - (30): 


A^= AXjj A^. -AZij A^.. 

Aq = cos £, (AZjj AAYjj-AYjj AAZjj ) + 

+ sin ^ (AXjj AAYjj-AYjj AAXjj ) 

a k = w d t sin E (AXjj A AZjj - AZjj AAXjj ) + 

+ sin % cos r| (AYjj A A Zjj - AZjj A A Yjj ) + 
+ cos E, cos T) (AXjj A AYjj - AYjj AAXjj )] 

Aoq = ~A K /W c j 

The equations (34) to (36) with 

sin £, s sin q = 0 
cos E, = cos T) = 1 


yield in matrix form 


A^ 


'-AZjj 0 AXj,] 

At) 

. AK/W d _ 


ff 


A AXy 

AAYy 

AAZjj 


(34) 


(35) 


(36) 

(37) 


(38) 


( 39 ) 
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Since there exist linear dependencies among the partial derivatives, it is not possible 
to estimate all of the parameters of interest from ground-to-ground baseline VLBI time 
delay observations. Of the initial 10 parameters of interest only six may be estimated 
simultaneously. 

It is obvious that the system of normal equations N = A^PA of the least squares 
adjustment will be singular. Two reasons have to be mentioned for this: 1) rank 

deficiencies due to coordinate system definition, and 2) rank deficiencies due to a critical 
configuration. 

There are certain conventions to overcome rank deficiencies of a ground-based 
geodetic VLBI network due to coordinate system definitions, see Bock (1980), Brouwer 
(1985), Dermanis and Mueller (1978), and Papo and Saleh (1988). Since in a Euclidean 
space, the coordinate system definition requires seven parameters; therefore, in practice the 
most common way to define the coordinate system is to constrain seven well-chosen 
“conventions” in the least squares adjustment. One possible choice for these seven 
parameters for ground-based geodetic VLBI is [Brouwer, 1985]: 

• the X,Y and Z coordinate of one station to define the origin (three translations) 

• the epoch ephemeris pole position to define the equatorial plane (two rotations) 

• the right ascension of one source to define the orientation in the equatorial plane 
(one rotation) 

• the velocity of light in vacuo as a scale parameter 

An other convention for the coordinate system definition may be the use of a 
minimum norm solution. This procedure yields a minimal trace for the variance-covariance 
matrix so that smoother looking variances appear (for details see Chapter 4.1). Here we 
mention that there exists a rectangular matrix G of m rows and r linear independent 
columns such that AG"^ = 0 . The result is that a nullspace L of N = A^PA is described. 
The basis of the nullspace is formed by the column vectors of the G matrix, the dimension 
of the nullspace is r . Such a matrix for the ground-based geodetic VLBI is given by, e.g., 
Brouwer (1985; Appendix E) and Dermanis and Mueller (1978). 

In order to circumvent the estimability problem, the earth rotation parameters (ERP: 
T), k) as one of the possibilities and right ascension (oq) may be redefined by introducing 
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A^ 0 = ^o 

(40a) 

Ario = ri -Bo 

(40b) 

Akq = k - Kq 

(40c) 

Aa 0 , = o-i - a 0 , 

(41) 


where £ 0 , Bo, k 0 are the adopted ERP providing the initial reference orientation of the 
network (CTS) with respect to the true equator and equinox, and a G is the fixed right 
ascension of one radio source (3C273B or (3 Persei) providing the reference origin of right 
ascension. Therefore, the corresponding estimable parameters are the ERP differences 
(A^o. Ario* Ak 0 ) as changes in the network orientation relative to the initial orientation and 
the right ascension differences Aa 0 ^with respect to the fixed right ascension a 0 . However, 
the errors in the four adopted basic parameters of orientation bias the baseline components. 
Thus, from this point of view, the baseline components AXjj , AYjj , AZjj are nonestimable 
and are replaced by the corresponding set of estimable components AX'jj , AY'jj , AZ'jj , 
contaminated by the above errors and defined by the following equations [Arnold, 1974], 
[Bock, 1980] and [Moritz and Mueller, 1987] 

"dAXyl r dAX'ij 1 f-AZij 0 AYjj IT d^ " 

dAYjj = dAY'jj - 0 AZjj -AXjj dTi 0 , 

dAZjj dAZ'jj AXjj -AYjj 0 dp G 

J J 1J J J (42) 

where dAXjj , dAYjj , dAZjj are the differential changes of the baseline components and 

dp 0 = doto - Wd dKo, (43) 

implying that the two differential rotations doto and dKo are inseparable, 
d^ 0 , dr] 0 , and dp o are the small errors of the initial reference orientations. Note that the 
baseline length is unaffected by the errors in the reference orientation being invariant of 
coordinate system definition. 

A list of estimable parameters recoverable from delay observations is: 

(AX jj, AY jj, AZ jj, 8;, Aa 0 ^, A^o, At|o> Ak g , ACq-, AC]-} (44) 
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A number of limiting conditions should be considered in the geodetic analysis phase 
of VLBI observations, among others some critical configurations. These are the result of 
an insufficient measurement design that yield a (nearly) singular system of normal 
equations in the least squares adjustment. The following two equations are useful to 
describe critical configuration situations. 

Rearranging Eq. (39) we arrive at the equation 


AZjj A k = W d (AYjj A^ + AXjj A n ) . 


(45) 


Furthermore, under the assumptions in Eq. (38) , 

Afy = AZjj AaZjj ctg5/ - 

- ( AXij + AY ij A A Yij ) *g (46) 

can be derived which with Eq. (45) is appropriate for studying critical baseline 
configurations, cf., [Bock, 1980], [Dermanis, 1980] and [Moritz and Mueller, 1987], A 
well-known example of critical configuration in ground-based geodetic VLBI that yields a 
singular case is the single baseline experiment with observed sources at only one 


declination, see Eq. (46). 

Recognizing that 

cos £ = cos q = 1, (47a) 

sin £> = %, (47b) 

sin T[ s T[, (47c) 

and neglecting products ^T|, an examination of Eqs. (1) and (20) reveals that since 

dijk* = Ki cos (0k - a t ) + K 2 sin (0k - a t ) + K 3 + K 4 (tk - to) (48) 

where 

Ki = - (AXjj - £, AZjj) cos 8 / (49a) 

K2 = (AYjj + r\ AZjj) cos 8/ (49b) 

K3 = -(AZjj + £, AXjj - q AYjj) sin 8/ + c ACo- ( 49c ) 

K4 = c ACijj (49d) 
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or more simply 

dijk/ = K sin (co e t k +<()) + K 3 + K4(tk - to) (50) 

where 

9 9 1/2 

K = (Kf + K|) (51a) 

<t> = a/-ao (51b) 

oco is the right ascension of the baseline at the initial epoch to, 

a time delay observation of a single ground-to-ground baseline using a single source 
represents a sinusoid superimposed on a straight line. Eq. (48) and (50) may help us in 
obtaining a better understanding of the dependencies among the unknown parameters. 

Counselmann and Shapiro (1978), Bock (1980), Dermanis and Grafarend (1981), 
Shapiro (1978), and Moritz and Mueller (1987) show which parameters are estimable from 
ground-based VLBI measurements. Of the geodetic parameters only baseline lengths and 
source declinations are estimable from time delay. All other parameters are variational in 
nature, that is they are determined relative to the initial orientation of the inertial and 
terrestrial reference frames. Only their variations (changes) are estimable. Therefore the 
initial values for right ascension and earth rotation parameters (polar motion components 
and UT1-UTC value) must be held fixed at a priori values with the use of external 
information (e.g., derived from an independent technique). Geocentric site coordinates are 
not estimable because the observations are primarily a function of the differences of these 
coordinates for each baseline. In addition these baseline components are contaminated by 
errors in the reference orientation, polar motion and UT1-UTC variations. 


Note that from one baseline time delay observations, only two of the earth rotation 
parameters are estimable. Therefore, multi-baseline configurations are needed to estimate 
all three earth orientation parameters. 


3.4.2 Time Delay Rate Model. Differentiating Eq. (20) with respect to time 

we get 


dijk* = co e {AXy [ cos cos 5* sin (0k - 0 ^) + sin ^sin r\ cos 5* cos (0k - a/)] 
+ AY;j cos r\ cos 8 1 cos (0k - a /) - 

- AZjj [sin cos 8* sin (0k - a i) - cos \ sin T| cos 8/ cos (0k - cx^)]} + 


+ c ACijj , 


(52) 
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where 


0) 



(53) 


is the spin rate of the earth, co e the instantaneous earth rotation vector. 

One sees that time delay rate is effectively insensitive to the AZjj component of the 
baseline, consequently, only the equatorial projection of the baseline can be estimated. In 
addition, the delay rate is uneffected by clock offset variations, AC oi j. 


Furthermore, since 

f = “exB (54 ) 

is orthogonal to co e and, thus, the origin of declination is undefined as well as the right 
ascension origin. 

Taking the differential of d ijW with respect to the parameters 
d(dijk/) = Bp dP p , 

p 

where the index p stands for the unknowns 

{AXy, AYjj, AZy, a t , 8 /, T\, k , AQj} 

The Bp’s in Eq. (55) are the required partial derivatives of the time delay with 
respect to the parameters subscripted p as follows: 


(55) 

(56) 


BAXij = <0 e [cos £, cos 8 f sin (0 k - o^) + sin % sin T| cos 8* cos (0 k - a ( ) ( 57 ) 

BAYij = w e COS T] COS 8* COS (0 k - CLf) (58) 

BAZjj = — co e [si n % cos sin (0 k — ct^) — cos £, sin q cos 8 1 cos (0 k — ot;)] (59) 

= co e {AX;j [ cos ^ s i° 'H cos 8; cos (0 k - a i) - sin % cos 8* sin (0 k - a^)] - (60) 

— AZ;j [cos £, cos 8^ sin (0 k — <Xi) + sin £, sin q cos 8/ cos (0 k - cq)]} 
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( 61 ) 


= co e [AXy sin E ) cos r\ cos 8/ cos (0k - a*) - 
- AYjj sin r( cos 8^ cos (0k - a^) + 

+ AZy cos E, cos T) cos 5/ (0k - a^) ] 


B k = co e Wd {AXy [cos E, cos 8 1 cos (0k - at) - sin ^ sin rj cos 8/ sin (0k - a/)] - 

- AY;j cos rj cos 8* sin (0k - a /) - (62) 

- AZjj [sin % cos 8^ cos (0k - a^) + cos E, sin r| cos 8/ sin (0 k - a^)]} 


B a ^ = WefAXjj [sin E, sin r| cos 8^ sin (0k - a^) - cos E, cos 8/ cos (0k - a*)] + 
+ AYij cos r| cos 8/ sin (0k - ct/) + 

+ AZy [sin E, cos 8; cos (0k - cq>) + cos E, sin ti cos 8/ sin (0 k - a/)]} 


(63) 


B§£ = - co e sin 8f {AXy [ cos £, sin (0k - a*) + sin E s sin T| cos (0k - a /)] + 
+ AYjj cos r| cos (0k - O-i) - 

- AZy [sin E, sin (0k - a^) - cos % sin rj cos (0k - a/)]} 
b ac 1;j = C . 


The following linear relationships are evident 
B^ = AXjj B AZjj - AZjj B AXi j 

B-q = COS % (AZjj B AYij - AYij B^y) 

+ sin E, (AXjj B AYij - AYij B AYy ) 


B K = W d 


sin rj (AXjj B AZy - AZjj B AXij ) 


+ sin % cos ri (AYjj B AZy - AZjj B AYjj ) 



B a/ = - V Wd 


(64) 

(65) 

( 66 ) 

(67) 

( 68 ) 
(69) 


indicating that it is impossible to estimate all of the parameters of interest from ground-to 
ground baseline VLB I delay rate observations. Eqs. (66) - (69) are formally very similar to 
Eqs. (34) - (37). 

The Eqs. (66) - (68) with assumptions in Eq. (38) yield in matrix form 
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( 70 ) 


B s 



BAXjj 



-AZjj 0 AXjj 


Bn 

= 

0 AZjj -AYjj 

B AYjj 

i 

CO 

JP 

1 


O 

5 

* 

*5 

1 

B AZjj 


The discussion of the parameters estimable from delay rate is identical to that of 
time delays except that in this case AZjj (AZjj) and ACq- are deleted, and declination 
differences 

A5 0 , = 5,-5o (71) 

replace 8/ . Thus, the list of the estimable parameters recoverable from delay rate 
observations is 

r 1 (72) 

■% AXjj, AY;j, AS 0 *, Aa 0 *, A^o> At|o> AKq , ACijj j • 


It can easily be shown that the equatorial baseline components are also contaminated by the 
errors of the initial orientation in the following sense 


dAXjj 


dAXjj 

r 

[dAYjjj 


dAYy' 



-AZjj 

0 


0 AYjj 
AZjj -AXjj 


d^ 0 

dri 0 

dp 0 


(73) 


where d(3o equals with Eq. (43) implying again that the two differential rotations doto and 
dKo are inseparable. 

Similar derivation to the Eqs. (45) and (46) gives 

AZjj B k = W d (AYjj + AXjj B n ) (74) 

B 5 , = - tgS; (AXjj B A Xij + AYjjBAYjj + AZjj B^jj) (75) 


Eqs. (74) and (75) can be used to detect critical configurations and for sensitivity analysis 
as well. 
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An examination of Eq. (52) or differentiating the Eq. (50) with respect to time, we 
get 

^ijW = R C0 e COS (CO e tk + ({)) + K 4 , (76) 

which is useful for an analysis. The coefficients K and K 4 are the same expressed by Eqs. 
(51a) and (49d). 

The delay rate model includes a reduced parameter set in comparison with the time 
delay model. The third component of the baseline is non-estimable and only declination 
differences may be estimated, see [Bock, 1980] and [Moritz and Mueller, 1 987]. 
Therefore, it is not possible with delay rates alone to estimate all of the parameters of 
geodetic and astrometric interest. 


3.5 Mathematical Models for Ground-to-Space VLBI Observables 

We are interested in forming equations which relate to observables to the unknown 
parameters suitable for least squares covariance analysis. Our main purpose is not to derive 
explicit observation equations, but rather to develop qualitative expressions to demonstrate 
the relationship between the observables (time delay and delay rate) and the parameters and 
to explore estimability problems inherent in the space VLBI system. 

3.5.1 Time Delay Model. From Eq. (9) we can get a simplified expression 
(similar to Eq. (16) [Pavlis, 1986]) suitable for the analysis as follows (the T is the 
transpose sign), see Fig. 13: 



X; 


Zj. 


R 2 (-5) R 1 H 1 ) R3(0k) - 


Yl 

4J 


IT 


+ c ACo r j + ACf rj (t k - to) 


cos 8 / 
sin 8 ; 
sin 8 / 


cos a 1 
sin a/ 


+ 


(77) 
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orbit 


Earth 


Fig. 13 Space VLBI time delay observable 


where 

Xj, Yj, Zj are the Earth-fixed coordinates of the station Pj , 

X[, Y£, Z{ are the coordinates of satellite S 1 in a true-of-date geocentric inertial 

coordinate frame at the epoch t k , 

a i, 5 1 is the right ascension and declination of the radio source in the same as 

above true-of-date coordinate frame, 

£ 5 r| are the polar motion components that relate the instantaneous rotation axis of 

the Earth with the average terrestrial pole, 

0 k is the Greenwich Apparent Sidereal Time (GAST) at epoch tk, 

c is the speed of light, 

to is the initial epoch of observation, 


ACo r j , ACi rj is the clock offset and drift between the reference clock at the 
telemetry/control station and the clock of station Pj, and 

Ri(<J>) is the rotation matrix for a right handed rotation <J> about the axis i . 
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The expression of Eq. (77) may be rewritten as 

djW = - Xj [cos £, cos 8/ cos (0k - (x^) - sin £, sin ri cos 8/ sin (0 k - a*) + 

+ sin % cos T) sin 8*] + 

+ Yj [cos T] cos 8; sin (0k - a*) + sin ri sin 8*] + 

+ Zj [sin £, cos 8^ cos (0k - ct-t) + cos \ sin rj cos 8; sin (0k - o^) - 

- cos £, cos rj sin 8^] + 

+Xj[ cos 8/ cos a* + 


+ Y][ cos 8/ sin a i + 

+ Zj[ sin 8/ + 

+ c ACo r j + AC{ r j (tk - to) 


( 78 ) 


In the following, the Cartesian coordinates Xj[ , Y£, of satellite S 1 are 
transformed into Keplerian orbital elements consisting of the semimajor axis a, the 
eccentricity e specifying the elongation of the orbital conic section, the inclination i 
specifying the orientation of the satellite’s orbital plane with respect to the equator of the 
Earth, the right ascension £2 of the ascending node, i.e., the angle measured eastward 
along the equator between the vernal equinox and the point where the satellite crosses the 
equator traveling in a northerly direction, the argument to of perigee, i.e., angle between 
the ascending node and the perifocal point measured positive with increasing mean 
anomaly, and the mean anomaly M , i.e., the sum of the mean anomaly at epoch and the 
product of the mean motion and the elapsed time from epoch. A transformation of Kepler 
elements 


{a, e i, co, £2, M} 


(79) 


into geocentric orthogonal coordinates {X,Y,Z} is given by the following formulas 
(dropping the lower and upper indexes) [El’yasberg, 1967; p. 63] and [Grafarend and 
Livieratos, 1978]: 


X 


cos £2 cos u - sin £2 sin u cos i 

Y 

= r 

sin £2 cos u + cos £2 sin u cos i 

_Z_ 


sin u sin i 
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where 


r = a (1 -e cos E) 

(81a) 

u = CO + f 

(81b) 

tg f/2 = Y(1 + e)/(l-e) tg E/2 

(81c) 

M = E - e sin E 

(8 Id) 

sin f = V 1 - e 2 sin E/(l - e cos E) 

(81e) 

cos f = (cos E - e)/(l - e cos E) 
E is the eccentric anomaly. 

(8 If) 


Equation (78) with Eq. (80) [and Eq. (77)] expresses the functional relationship of 
the time delay observation with the listed parameters. Of direct geodetic interest are the 
station coordinates Xj, Yj, Zj , the Keplerian orbital elements Eq. (79) of the orbiting 
radiotelescope S 1 and the earth rotation parameters (ERP) 5, T|, k . The source coordinates 
(xi and 8/ are of astrometric interest, while the clock offset ACo r j and rates ACi r j are 
nuisance parameters defined to make the mathematical model more realistic. 

The GAST, 0 k as well as the Earth rotation parameters K are reformulated 
identically to that of ground-to-ground baseline VLBI observables, see Chapter 3.4, i.e., 
Eqs. (17) and (23). Assume that for an observation session (one day or one revolution of 
Keplerian orbit) only one set of earth rotation parameters (^, Tj, k) is determined. Thus, 
our parameter set will contain the following 16 parameters: 

(Xj, Yj, Zj, a, e i, (0, Q, M, £,, T|, K, a e , 6 1 , ACo 1 ^, AC{ rj } (82) 

In the following we examine the complete observation eqation Eq. (78) with Eq. 
(80) for the parameters, the coordinates of the ground VLBI station Pj, the Keplerian orbital 
elements of the orbiting radiotelescope S', the ^th radio source positions, the earth rotation 
parameters, and the clock parameters to determine which of them are separable or 
nonseparable and, therefore, design experiments in which those parameters of primary 
interest will be estimable. 

Since Eq. (78) is non-linear with respect to that parameter set, we proceed as usual 
to linearize Eq. (78). Taking the differential of dj w with respect to the parameters, 

i(^) = ICpdP p . (83) 
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where the C p ’s are the required partial derivatives of the time delay with respect to the 
parameters of interest, indexed by p. They are as follows: 


C X j = - [cos \ cos 5^ cos (0 k - a/) - sin % sin r\ cos 8f sin (0 k - a/)] 

+ sin % cos T) sin 5;] 


(84) 


Cyj = cos r\ cos sin (0 k - a/) + sin rj sin 8^ 

Czj = sin % cos 8/ cos (0 k - a^) + cos % sin rj cos 5* sin (0 k - a*) - 

- cos % cos T) sin 5/ 


(85) 

( 86 ) 


= Xj [sin t, cos 8/ cos (0 k - a i) + cos £, sin r\ cos 8* sin (0 k - a /) - 
— cos % cos T| sin 8/] + 

+ Zj [cos E, cos 8* cos (0 k - ai) - sin £, sin r| cos 8/ sin (0 k - a/) + 

+ sin E, cos rj sin 8^] 

Ct| = Xj [sin E, cos T) cos 8^ sin (0 k - Ot*) + sin E, sin TJ sin 8/] - 
- Yj [sin r\ cos 8^ sin (0 k - a^) - cos r\ sin 8*] + 

+ Zj [cos ^ cos r| cos 8 1 sin (0 k - a/) + cos E, sin r| sin 8*] 


(87) 


( 88 ) 


C K = Wd {Xj [ cos E, cos 8* sin (0 k - a^) + sin E, sin rj cos 8^ sin (0 k - a/)] + 

+ Yj cos T) cos 8^ sin (0 k - a i) - 

-Zj[ sin E, cos 8i sin (0 k - a/) - cos E, sin T) cos 8/ cos (0 k - o^)]} 

C a = (1 - e cos E) {cos 8^ [cos u cos (Q - ai) - 
- sin u cos i sin (£2 - a/) ] + sin 8^ sin u sin i} 


(89) 

(90) 


C; = r [sin u sin i sin (Q - a^) cos 8^ + sin u cos i sin 8*] 

C w = - r {cos 8/ [sin u cos (Q - at) + cos u cos i sin (Q - 8^)] 

- sin 8^ cos u sin i} 

Cq = - r cos 8/ [cos u sin (Q - 8^) + sin u cos i cos (£2 - 8;)] 


a(e-cosE) c + 
(1 - e cos E) 2 


sin f 


+ 


sin f 


1 - e cos E 1 - e 2 


p. , _ a e sin E p , sin f P 

'-M — ~ '-a "l : ^CO 

(1-ecosE) 2 sinE (1 -ecosE) 


(91) 

(92) 

(93) 

(94) 

(95) 
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( 96 ) 


Co = - {Xj [cos t, cos 8* sin (9k - «/) + sint, sin ri cos be cos (9k - a/)] + 
+ Yj cos r| cos 8/ cos (9k - ae) - 

- Zj [sin £ cos 8/ sin (9k - a*) - cos£ sin rj cos 8* cos (9k - a*)]} + 

+ r cos 8* [cos u sin (£2 - a e ) + sin u cos i cos (£2 - a/)] 

C 5i = sin 8 1 {Xj [cos £ cos (9 k - a t ) - sin£ sin t\ sin (9 k - a/)] - 
- Yj cos rj sin (9 k - a/) - 

- Zj [sin % cos (9k - a /) + cos£, sin T| sin (9 k - a*)] - 

- r [cos u cos (Q - a e) - sin u cos i sin (£2 - a*) ]} - 
be [Xj sin ^ cos Tj — Yj sin rj + Zj cos % cos r| - r sin u sin i] 


cos 


CACo-j - C 

CaC, i . = C (t k — to) 
hj 


(97) 


(98) 

(99) 


In order to get the partial derivatives of time delay (dj^ = d) with respect to 
Keplerian orbital elements, the following relationship: 


C a = 


Q = 


3d 

3a 

3d 

3e 


r - — 

{ ~ i ~ d i 

r _*L 

'-CO - 

3co 

3d 

dQ 

3d_ 
3M 


Ca = 


Cm = 


3X 

3Y 

3Z 

3a 

3a 

3a 

3X 

3Y 

3Z 

3e 

3e 

3e 

3X 

3Y 

3Z 

3i 

3i 

3i 

3X 

3Y 

3Z 

3(0 

3co 

3co 

3X 

3Y 

3Z 

3£2 

3£2 

3Q 

3X 

3Y 

3Z 

3M 

3M 

3M 


3d 

3X 1 

3d 

3Y 1 

3d 

3Z 1 


= Cvl 


— Cvl 


— C7I 


( 100 ) 


was applied, where the partial derivatives of time delay with respect to the orbiting 
radiotelescope’s cartesian coordinates are as follows: 


C x i = cos be cos ae 
C Y i = cos be sin ae 
C z i = sin be 


( 101 ) 

( 102 ) 

(103) 
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The elements of the coefficient matrix in Eq. (100), called as a linear Kepler transformation 
matrix K by Grafarend and Livieratos (1978), are given in Appendix A. We call this 
matrix here as a linear Kepler transformation matrix Kc of position (or coordinate) type. 


A useful check of the partial derivatives of the Cartesian inertial coordinates with 
respect to Keplerian orbital elements f 3x /a Kj l is their use with the corresponding partial 


derivatives of Keplerian elements with respect to the Cartesian parameters 
confirm the relation 

V 3XlM = 5ik , 


M 


dX k Jto 


j r 1 3Kj ax k 


(104) 


for all i, k e {1,2,3} 

where 8ik is the Kronecker-delta, X represents the Cartesian coordinates {X, Y, Z} and 
K represents the Keplerian orbital elements {a, e, i, co, Q, M) . Both types of partial 
derivatives in Eq. (104) are collected in Appendix A. To compute the partial derivatives of 
Keplerian elements with respect to the Cartesian parameters, the partials given in Hill 
(1989) are used. Appendix A also includes a computer program with subroutines as well 
as a practical example to confirm the relation in Eq. (104). 

The partial derivatives [Eqs. (84) - (99) and Eqs. (101) - (103)] of the time delay 
observable with respect to the corresponding parameters of the mathematical model Eq. 
(78) with Eq. (80) constitute the elements of the design matrix C in Eq. (83) that forms 
the normal matrix N = C^PC (P is the weight matrix of the observables). The inverse of 
that normal matrix yields the a priori covariance matrix of parameters. The magnitude of a 
particular partial derivative which reflects the sensitivity of an observation to a particular 
parameter determines its numerical contribution to the normal matrix. Most of the partial 
derivatives listed in Eq. (84) - (99) are diurnal sinusoids, some others are composed of 
diurnal sinusoids, functions of station coordinates, and Keplerian orbital elements (or 
varying satellite positions). 

Examining the analytical expressions for the partial derivatives we can derive some 
interesting sensitivity relations between the time delay observable and the parameters to be 
solved for. Since the partial derivatives of delay with respect to the polar components of 
the ground station and satellite position and to the declination of the sources (when Zj = Zl) 
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go to zero as the declination goes to zero, the sensitivity of delay to these parameters 
decreases as the declinations go to zero. The tendency of the sensitivity to the right 
ascension of the sources and the earth rotation parameters K is just opposite, that is, from 
Ck and C u ^ it is evident that the time delay is sensitive to the declinations of sources near 
the equator. The orbiting radiotelescope’s coordinates (Keplerian elements) effect only the 
partial derivatives of delay with respect to radio source position components, see Eqs. (96) 

- (97). From the analytical expression of these partial derivatives (C a( and C~ bf ) we can 

gather that their numerical values will in general increase as the coordinate separation 
between the station and the orbiting radiotelescope increases (advantage of space VLBI). 

With three coordinates for the ground station, six Keplerian orbital elements for the 
orbiting radiotelescope, two parameters for the radio source position, three for the ERP and 
two for the clock, there are 16 parameters in total in the model Eq. (78). Therefore, the 
design matrix C given in Eq. (83) will have a minimum of 16 columns, if one ground 
station and one radio source is considered. In identifying linear dependencies between the 
parameters, we can separate C into five submatrices according to the five major groups of 
parameters, namely, the station coordinates, the orbital elements, earth rotation parameters, 
radio source positions and clock parameters. From Eqs. (84) - (99) we can state that there 
are no linear dependencies between columns of the design matrix within the same group. 
However, dependencies do exist between columns of the design submatrices among 
different groups with the exception of the ones belonging to the clock parameters. The 
columns of the station coordinates can be written as linear combinations of those of the 
earth rotation parameters. Therefore, we cannot separate the ground station parameters and 
the ones for the earth rotation in a simultaneous adjustment. 

Furthermore, one can easily write the linear relationship between columns of the 
station coordinates, the right ascension of radio sources and the longitude of satellite orbit’s 
ascending node (Q ) . The linear dependence between the columns of the right ascension 
(a) and earth rotation parameter (k) is also obvious. Consequently, the following linear 
relationships are evident among the partial derivatives: 

Q^XjCzj-ZjCxj (105) 

C n = cos^(Z j C Yj -YjCz j ) 

+ sin B, (Xj C Yj - Yj C Xj ) ( 106) 
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C K = W d [simi(XjCz j -ZjCx j ) 

+ sin % cos r| (Yj C Zj - Zj C Yj ) 

+ cos E, cos r| (Xj C Y j - Yj C Xj ) 

= - °Vw d - 

Recognizing that sin ^ = sin T| = 0 and cos E, = cos T| = 1 , the linear relationships through 
Eqs. (105) and (107) are reduced to the following matrix form which is more clearly 
indicating the dependencies among the corresponding partial derivatives: 

r 

Cn 

c ‘/w d 

The linear dependencies derived analytically by Eqs. from (105) to (108) are 
formally very similar to ones existing among the partial derivatives of ground-based time 
delay observable (see Eqs. (34) - (37)). It is quite obvious from a comparison of Eqs. 
(109) and (39). The only difference is the use of station coordinates instead of baseline 
components. However, a comparison of the Eqs. (37) and (108) shows that in the latter 
equation a new term (Cq) appears due to the additional linear dependency of radio source’s 
right ascension with the ascending node £2 of the satellite’s Keplerian orbit in the space 
VLBI. A geometric interpretation of this reason is that all three parameters, i.e., right 
ascension a t , Greenwich Apparent Sidereal Time (0k) (or earth rotation parameter k ) 
and ascending node Q of the Keplerian orbit are “measured” in the same plane, that is, in 
the equator (see Fig. 14). 

In the following, it is useful to present the geometric interpretations of the rank 
deficiencies expressed in Eqs. (109) and (108). Eq. (109) shows a linear dependence 
between various combinations of ground station coordinates Xj, Yj, Zj and the earth 
rotation parameters q, k . These indicate a rank deficiency of three due to a lack of 
absolute orientation of the network of ground stations with respect to the true-of-date 
celestial frame which cannot be sensed by the observables. Eq. (108) is due to a lack of 
reference direction (origin or right ascension of the radio source and of the ascending node 



(107) 

(108) 
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of satellite orbit), for the true-of-date celestial frame, i.e., the observables are insensitive to 
the orientation of the true-of-date inertial frame in right ascension (see Fig. (14)). 





Fig. 14 A geometric configuration of the Greenwich Apparent Sidereal Time (GAST), 
the right ascension (a£) of a radio source and the ascending node ( Ci.) of the 
Keplerian orbit. 


It has been shown in Chapter 3.3 that the time delay observables (or time delay 
times the speed of light as the path difference) are invariant under datum translations and 
rotations, but not invariant under datum scale change. The observable is insensitive in the 
location and orientation of the coordinate system where the coordinates of the end points of 
the baseline are given. Therefore, time delay observations cannot provide datum origin and 
orientation, but only scale. The missing origin and orientation components are the datum 
defect. 


There are principally two reference frames associated with the mathematical model. 
The first is an earth-fixed frame, and the second is a true-of-date inertial frame. Since time 
delay can only provide scale, nine parameters define the two frames — three for the origin 
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source 


* 



Fig. 15 Geometry of time delay observable 

and three for the orientation of the earth-fixed frame, plus three for the orientation of the 
true-of-date inertial frame with respect to the earth-fixed frame. The design matrix C 
contains partial derivatives of time delay with respect to parameters which include station, 
satellite and radio source coordinates. Because such coordinates are not estimable 
quantities, the design matrix should have a column defect of at least nine, provided there 
are no configuration defects. 

However, the time delay observable will remain unchanged if the wavefront rotates 
in such a way that it remains tangent to a cone with axis the baseline and half angle equal to 
\|/ (Fig. 9 and 15), see [Counselman and Shapiro, 19781. From that fact, two additional 
degrees of freedom exist (three rotations of a plane minus one constraint of tangency to the 
cone). Furthermore, in practice, instead of three coordinates we need only two directional 


49 



coordinates of source, that is, the right ascension oc and declination 8 to define the true- 
of-date inertial reference frame. Thus, seven rank deficiencies should appear in the 
mathematical model Eq. (77) of the observable. 

We have shown that the three translational and the three rotational parameters 
cannot be found from the either ground-based and space VLBI time delay observations. 
There are certain conventions in ground-based VLBI to define the origin, that is, e.g., to 
fix the X,Y and Z coordinates of one station or to use a minimum norm solution. In the 
latter case, we will have an optimal datum definition in the classical sense with the help of 
the standard inner constraints for 3D networks [Brouwer, 1985; Appendix E], 

Note that we were searching to analytically find linear constraints for “translational” 
datum defect for space VLBI observation equations. In the case of ground-based VLBI, it 
is easy to get a nullspace basis in such a way that AG?" = 0 for translational defect in the 
time delay observation equation [Brouwer, 1985] and [Dermanis and Mueller, 1978]. 
However, a comparison of Eqs. (13a) and (13b) shows that in the case of space VLBI, the 
derivation of such a G? - matrix (CG? = 0) for definition of the origin in the minimum norm 
solution is impossible. Therefore, it implies that some fixed values should be applied to fix 
the origin of the coordinate system for space VLBI observables. 

Contrary to ground-based VLBI network, in the space-VLBI the origin of the 
reference frames is fixed by the given Keplerian orbital elements at some reference epoch, 
T 0 . In fact, the orbital elements are not time invariant, i.e., their derivatives with respect to 
time, Ei , will not be zero. The variation of a given element from some reference epoch, 
T 0 , to the epoch of utilization, T , can symbolically be described by the following equation 

Ei = E ‘ D + J ^ dT (HO) 

where EJ’ is the element in question at the initial epoch T 0 , and Ej at T (Fig. 16), 
see [Moritz and Mueller, 1987; pp. 368-369]. The integral represents the perturbation of 
the element Ej . 

A rank deficiency of a normal equation system may easily occur as a combined 
effect of coordinate system definition and critical configurations which are the result of an 
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insufficient measurement set-up. Therefore, the consideration of critical configurations is 
also very important. 

From Eq. (109) one obtains 

Z j C K = W d (Y j q + X j C n ) (Hi) 

Considering the assumptions in Eq. (38) 

Cs t = ctg b ( Zj C Zj - tg §t (Xj C X j + Yj C Yj + a C a ) + r sin u sin i/cos 8 * (] j 2) 

is derived. In searching of the rank deficiencies the following equation 



Fig. 16 Keplerian orbital elements: Right ascension of node, Ej ; inclination, E 2 ; 

argument of perigee, E 3 ; semidiameter of the elliptic orbit, E 4 ; eccentric of the 
ellipse, E 5 ; time of perigee passage, E 6 
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cos i C ro = ctg u sin i Q + C& 


(113) 


has been obtained. Eqs. (Ill) - (113) can be used for studying the critical station and 
satellite configurations. They may be useful in observation design works as well. For 
instance, additional rank defect occurs by a special choice of Keplerian orbital elements. If 
i = 0° , then from Equation (1 13) it follows that 


Cco = Q2 (114) 

indicating an inseparable situation for the parameters (0 and £2 . In the space VLBI 
network this special situation gives a rank defect of internal type [Graferend and Livieratos, 
1978] 


Because of the linear relationships Eqs. (108) and (109) among the coefficients Cp 
of the error Eq. (83), it is not possible to determine all of the 16 unknowns in Eq. (82). 
Therefore, a new set of estimable parameters must be introduced. Concerning the linear 
relationships expressed in Eqs. (105) - (107) [or see Eq. (109)] and in Eq. (108), the ERP 
differences Aq 0 , Apo , Ak 0 and right ascension difference Aa 0 / may be introduced as 
the corresponding new estimable parameters like that of the ground-based VLBI time delay 
observable, see Eqs. (40) - (41). However, the errors in the four adopted basic parameters 
of the initial orientation (£ 0 , r| 0 , Kq, a 0 /) biase the station coordinates. Thus, from this 

point of view, the station coordinates Xj U , Yj U , Zj U are non estimable. They are replaced 

/ / / 

by the new set of estimable station coordinates Xj, Yj, Zj contaminated by the errors of the 
initial orientation (d^ G , dri 0 , dKo, doc 0 ) . It can be shown that the errors of the 
corresponding parameters can be expressed in the following manner: 


dXj 1 


dXj 


' -Zj o Yj ' 


d^o 

JlSi 

i 


“S’ ¥ 

1 


0 Zj -Xj 
Xj -Yj 0 _ 


dr|o 
. dpo . 


where dXj, dYj, dZj are the differential changes of the station coordinates and dp Q is 
identical with Eq. (43), that is, the two differential relations doto and dK 0 are inseparable. 

Furthermore, the error of the initial reference orientation in right ascension makes 
its influence in the space coordinates of the orbiting radiotelescope. Consequently, the 
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satellite coordinates (only the equatorial ones, x£ and y£ ) will be contaminated by the 
error da*, in the following manner: 




r 

S 

1 


0 

a 

x) 

0 


Xt 

1 


1 


— da 0 0 


Y^ 


where dX£ and dY| are the differential changes of the orbiting radiotelescope’s equatorial 
coordinates. 


Note that the declination of the reference source should be nearly equatorial to 
provide a strong definition, similarly to that of the case of ground-based VLBI. This can 
be seen by an examination of C a/ , Eq. (96). 


In order to show the main differences between the ground-based and space VLBI 
time delay observables, under the assumptions of Eq. (47) we rewrite Eq. (78) in a similar 
way to that of Eq. (48) as follows: 

djw = Ki cos (0k - a t ) + K 2 sin (0k - a t) + K 3 (tk) + K 4 (tk - to) , qh\ 


where 


Ki = - (Xj - 1 , Zj) cos 6/ 

K 2 = (Yj + n Zj) cos 8; 

K 3 (t k ) = - (Zj + Xj £ - Yj Ti) sin 8/ + C ACl t} 

+ Xj cos 8 / cos a 1 
+ Y£ cos 8 1 cos a 1 
+ Z£ sin 8; 


K 4 = c ACi r j 


(118a) 

(118b) 

(118c) 


(1 18d) 


or more simply 

djk* = K sin (co e tk + (?) + K 3 (tk) + K 4 (tk — to) 


(119) 


where 0 is the phase of the sinusoids relative to some initial epoch, to e is the rotation 
rate of the earth, and 
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K - (K [ + K^) 1 ^ 2 


( 120 ) 


Eqs. (117) and (1 19) shows that a space VLBI time delay observation between a 
single ground station and the orbiting radio telescope using a single source represents a 
combined function of a sinusoid and a periodic function with a period of the orbital 

revolution. 

3.5.2 Time Delay Rate Model 

The geometric delay rate is the time derivative of the geometric delay. Including the 
clock parameters, the delay rate may be expressed by Eq. (10). Differentiating Eq. (78) 
with respect to time we get 

djw = w e {Xj [cos % cos 5/ sin (0 k - a*) + sin £, sin ri cos 6/ cos (0 k - a,)] + 

+ Yj cos r| cos 8 1 cos (0 k - a^) - 

-Zj [ sin £, cos 8/ sin (0 k - a /) - cos t, sin r| cos 8* cos (0 k - a*)]} 

+ X k cos 8* cos a t + 

+ Y k cos 8 1 sin QLt + 

+ zJc sin 8* + 

+ c AC} rj , (121) 


where 

co e is the spin rate of the earth expressed by Eq. (53) and 

X k Y k zJc are the velocity components of satellite in a true-of-date geocentric 
inertial coordinate frame at the epoch t k . 

Instead of the Cartesian velocity components X k , Y k , 2* of the orbiting radiotelescope S 1 , 
their equivalents expressed in Keplerian orbital elements are used by the following 
formulas (dropping the lower and upper indexes) [El yasberg, 1967, p. 63]. 


X = v r (cos Q cos u - sin Q sin u cos i) - 
- v u (cos Q sin u + sin Q cos u cos i) 

(122a) 

Y = v r (sin Q cos u + cos Q sin u cos i) - 
- v u (sin Q sin u - cos Q cos u cos i) 

(122b) 
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Z = v r sin u sin i + v u cos u sin i , 


(122c) 


where 

v r = ifjl/p esinf (123a) 

v u = Vm7p (1 +ecosf) (123b) 

p = a(l-e 2 ) (123c) 

jj. is the gravitational constant times the mass of the Earth . 

Our parameter set contains the same unknown parameters indicating by Eq. (82) 
with the exception of clock offset ACo rj to which the time delay rate is uneffected, 

{Xj, Yj, Zj, a, e, i, co, Q, M, rj, k, a t , 8<, ACi rj } q 24 ) 

Taking the differential of dL^ with respect to the parameters listed in the previous 

chapter 

d ( d lw) = X D P dP P - 

P (125) 

where the D p ’s are the required partial derivatives of the time delay with respect to the 
parameters subscripted p as follows: 

D X j = co e [cos E, cos 8/ sin (0k - a*) + sin E, sin r\ cos 8* cos (0k - a^)] (126) 

D Y j = <». cos rj cos 8/ cos (0k — o./) (127) 

Dzj = - co e [sin E, cos 8; sin (0k - a;) - cos E, sin T| cos 8 1 cos (0k -a;)] (128) 

= co e {Xj [cos £ sin r\ cos 8* cos (0k - a*) - sin i; cos 8/ sin (0k - a/)] - 

-Zj [cos E, cos 8; sin (0k - a;) + sin E, sin rj cos 8* cos (0k - a^)]} (P9) 

D p = co e [Xj sin E, cos t| cos 8^ cos (0k - a/) - 

- Yj sin rj cos 8; cos (0k - at) + 

+ Zj cos E , cos T) cos 8* cos(0k - a^)] (130) 
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(131) 


D k = co e W d {Xj [cos % cos 5/ cos (0 k - a t ) - sin % sin ri cos 8* sin (0 k - 8/)] - 

- Yj cos r\ cos hi sin (0 k - a*) - 

- Zj [sin % cos 8* cos (0 k - a :/) + cos % sin ri cos 8/ sin (0 k - 8*)]} 


D a = cos 8/ | - J [cos u cos (ft- a /) - sin u cos i sin (ft - a/)] + 

+ ^-[sin u cos (ft - a i) + cos u cos i sin (ft - a*)] f - 
2a J 

- sin 8/ Qk sin u sin i + ^- cos u sin ij 


Di = cos 8* [v r sin u sin i sin (ft - a*) + 
+ v u cos u sin i sin (ft - a*)] + 

+ sin 8/ ( v r sin u cos i + v u cos u cos i) 


D ro = - cos hi {v r [sin u cos (ft - a*) + cos u cos i sin (ft - a*)] + 
+ v u [cos u cos (ft - ai) - sin u cos i sin (ft - a*)]} + 

+ sin 8* ( v r cos u sin i - v u sin u sin i) 


(132) 


(133) 


(134) 


Djj = -cos 8/ {v r [cos u sin (ft - a*) + sin u cos i cos (ft - a/)] - 
- v u [sin u sin (ft - a /) - cos u cos i cos (ft - a/)]} 


(135) 


D e = cos 8/ f - v r — — * — [cos u cos (ft - at) - sin u cos i sin (ft - a/)] - 

e 1 e(l - ecosE) L 

_ Vll e ± cos f [ sin u cos (ft - a*) + cos u cos i sin (ft - a*)]} + 

(1 -e 2 ) (1 + e cos f) L 


+ V u 


+ sin 8/ 
e + cos f 

(1 -e 2 ) (1 +ecosf) 


Vr e(l- ecosE) 

} 


sin u sm i + 


cos u sin l 


(136) 


_ v r cos f v u sinf r cos r cos u cos (Q _ _ sin u cos i s j n (Q _ a< )] + 

sin E (1 -e cos E) 


sin 8; sin u sin i} . 


(137) 
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D a< = We {Xj [sin \ sin rj cos 8^ sin (9k - a^) - cos £, cos 8/ cos (0k - ct/)] + 
+ Yj cos rj cos 5/ sin (0k - a*) + 

+ Zj [sin £, cos 8t cos (0k - a^) + cos £, sin r) cos 5* sin (0k - a*)]} + 

+ cos 8; { v r [cos u sin (£2 - a^) + sin u cos i cos (£2 - a;)] - 
- v u [sin u sin (£2 - a^) - cos u cos i cos (£2 - a/)]} 

D§^ = - co e sin 5* {Xj [cos ^ sin (0k - a i) + sin £ sin T| cos (0k - a*)] + 

+ Yj cos r\ cos (0k - a/) - 

- Zj [sin ^ sin (0k - a /) - cos £, sin r| cos (0k - a*)]} - 

- sin 8i {v r [cos u cos (£2 - a i) - sin u cos i sin (£2 - a/)] - 

- v u [sin u cos (£2 - a^) + cos u cos i sin (£2 - a*)]} + 

+ v r sin u sin i cos 8/ + v u cos u sin i cos cos 8/ . 


( 138 ) 


(139) 


Dac! =C 
■In 


(140) 


In order to get the partial derivatives of time delay rate (djw = d) with respect to 
Keplerian orbital elements, the following chain-rule differentiation 


D = — 

a a a 


ax at az 


da 3a 3a 

ae 


ax at az 
ae ae ae 

0 1 


ax at az 


ai ai ai 

D - — 

dCO 


ax at az 
aco aco a® 

n - 3d 

Dn_ an 


ax at az 


a£2 a£2 a£2 

? 

II 

QJ 


ax at az 


3M 3M 3M 


a a 
ax 1 

ad 

at 1 


= D x I 
= DyI 


ad 

az 1 


Dz 1 


(141) 


was applied, where the partial derivatives of time delay rate with respect to the orbiting 
radiotelescope’s Cartesian velocity components are as follows: 

D x I = cos 8i cos OLi (142) 
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Dy 1 = cos 5; sin a i 
D £i = sin a i 


(143) 


(144) 


The elements of the coefficient matrix in Eq. (141) are derived and given in Appendix A. A 
useful check of the partial derivatives of the Cartesian velocity components with respect to 
Keplerian orbital elements 5X;/3Kj is their use with the corresponding partial derivatives 
of Keplerian elements with respect to the rectangular velocity components 8Kj/9X k to 
confirm the relation 



3*i^U 5lk 

3Kj dX k 


for all i, k e {1,2, 3} 


(145) 


where 8 ik is the Kronecker-delta, X represents the Cartesian velocity components 
{X, Y, Z} and K stands for the Keplerian orbital elements {a, e, i, co, Q, M) . 

The magnitude of the terms containing ^ and T\ in Eq. (121) are small, indicating that 
the delay rate is effectively insensitive to the Zj component of the ground station. More 
precisely, it is completely insensitive to the coordinate component of the ground station 
parallel to the Earth’s angular velocity vector. This is a disadvantage for geodetic 
applications of space VLBI delay rate observables. In addition, the delay rate is unaffected 

by clock offset variations, AC^ . From D^i it is evident that the delay rate is insensitive 

to the polar component of the orbiting radiotelescope’s velocity as the declination goes to 
zero. Only the partial derivatives of delay rate with respect to radio source position 
components contain the velocity components of orbiting radiotelescope. 

In the ground-based VLBI, the delay rate observable is insensitive to the declinations of 
sources near equator due to the fact that B§^ in Eq. (64) contains only the polar component 
of the declination. However, in the space VLBI the delay rate observable is more sensitive 
to the declinations of sources near the equator, since the corresponding partial derivative 
Dg / in Eq. (139) contains the equatorial component of the radio source’s declination as 
well. Therefore, an advantage of space VLBI over ground-based VLBI is that the space 
VLBI may be used for accurate determination of the declination of those radio sources near 
equator by the use of such ground-to-space baselines that possess large polar components 
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and relatively large values for the polar component of the radiotelescope’s velocity, 
compare Eqs. (64) and (139). 

The following linear relationships can be detected among the partial derivatives D p : 


J5? 

ii 

0 

dS 

1 

_N 

J? 

(146) 

Dt, = cos S, (Zj D Yj - Yj D^) + 

+ sin ^ (Xj D Yj - Yj D Xj ) 

(147) 

D k = W d [sin TJ (Xj D Zj - Zj D Xj ) + 

+ sin £, cos r\ (Yj Dzj - Zj D Yj ) + 
+ cos £, cos r| (Xj D Y j - Yj D X j) 

(148) 


D a ^ = - D k /W d - Dq (149) 

Eqs. (146) - (148) with assumptions in Eq. (38) yield in matrix form 






D X] 1 



-Zj 0 Xj 


J 

Dn 

= 

0 Zj -Yj 



D^W d J 


-Yj Xj 0 


PZj 


The relations from (146) to (150) are formally similar to the ones derived for the time 
delay model, see Eqs. (105) - (109). There is a formal similarity even with the analytically 
derived equations of the linear dependencies among parameters of the ground-based VLBI 
time dely and delay rate observation models. Compare the Eqs. (34) - (39), Eqs. (66) - 
(70), Eqs. (105) - (109) and Eqs. (146) - (150). 

From (150) we can write 

Z j D K = W d (Y J D^ + X J D 11 ) (151) 

Similar to Eqs. (112) and (113) for time delay observable, we can derive 

= - tg 8/ (Xj Dxj + Yj Dyj + Zj Dzj - 2aD a ) 

+ (v r sin u + v u cos u) sin i/cos 5^ (152) 
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and 


cos i D w = Da + sin i 


v r cos u - v u sin u 
v r sin u + v u cos u 


Di 


(153) 


Eqs. (151), (152) and (153) are useful to study the critical configurations. 


The discussion of the parameters estimable from delay rate is identical to that of time 
delays except that in this case Z) (Zj) and ACo rj are deleted. Initial reference orientation 
parameters "Ho, K o and Oq may be introduced again in order to overcome the datum 
rank defect problem. However, their errors effect the equatorial coordinates of the ground 
station and the equatorial velocity components of the orbiting radiotelescope in the 
following manner: 


dXj 

dYj 


dXj 

dY j 


-Zj 

0 


0 Yj 

Zj -Xj 


l[ d ^ 0 

drio 

J L dp 0 


(154a) 


where dp o is identical with Eq. (43), that is, the two differential rotations doco and diCo 
are inseparable and 


dX k 


dX k 


0 da G 


Xk 

dY k 


dY k 


-da 0 0 


Yk 


(154b) 


Differentiating Eq. (117) we get 

djk/ = K co e cos (co e t k + <t>) + K 3 Ok) + , ( 155 ) 


where 

K 3 (t k ) = X{ cos 8 1 cos ai 
+ Y k cos 8/ sin cti 
+ zl sin 8^ , 


(156) 
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showing that a space VLBI time delay rate observable represents a combined function of a 
sinusoid and a periodic function of satellite’s velocity components with a period of the 
orbital revolution. 


3.6 Mathematical Models for Space-to-Space VLBI Observations 


3.6.1 Time Delay Model. Eq. (11) may be rewritten as 


'Xl-xf 

cos 5/ cosa/ 


r 

1 

>- N 
l l 

1 

cos 8^ sina^ 
sin 8t 

+ c 



(157) 


where x£, Y^, z { andXj[,Yj[, Zj[ are the coordinates of the orbiting radiotelescope S 1 
and S J , respectively, in the true-of-date inertial reference frame at the epoch tk. 


In the model Eq. (157) we have ten unknown parameters 


|xL yL zl xl Yl Zl 8/, AC(jJ q , ACP rq J 


(158) 


to be determined. Taking the differential of dj^ with respect to the parameters, 

■fi-SEpdPp. 

p (159) 

where the E p ’s are the required partial derivatives of the time delay with respect to the 
parameters of interest, indexed by p. They are as follows: 


E x i = -cosS^ cosa; 

(160) 

Eyi = -cosS/ since/ 

(161) 

E^ = -sinS/ 

(162) 

E X J = cos8/ cosa / 

(163) 

E y J = cos8/ sinot/ 

(164) 

EzJ = sinS/ 

(165) 


61 



E a/ - cos 5* [(Xjf - X{) sina t - (Yj - Yj) cosa* 

E§^ = sin5/ [(Xfc - X£) cosa t + (Yt - Yi [)] since* - (Z{ - Z{) 
E*# =c 

u rq 

E acF = c ( t k“to) 

irq 

The following relationships are immediately apparent from above set 

ExJ^-Exj' 

E Y j‘=- E Y i ' 

B 4 =_Ez k 


These parameters cannot be separated at all, and new combined unknowns must be 
introduced: 


> 

X 

II 

X 

1 

X 

(173) 

ay£=y£-y£ 

(174) 

> 

is, 

11 

£ 

1 

JKL 

(175) 

with partial derivatives 


E<=E X I 

(176) 

E A y” = Ey k 

(177) 

1 

II 

(178) 


(166) 

c °s5/ (167) 

(168) 

(169) 

(170) 

(171) 

(172) 


The new parameters are coordinate differences of the orbiting radiotelescopes S 1 and S J at 
the epoch tk . 

If the Keplerian orbital elements are used instead of the satellite’s Cartesian 
coordinates, then the partials in Eqs. (90) - (95) may be used for the formulation. In this 
case, it can be shown that the ascending nodes (Qi and Q 2 ) of these two satellites (S 1 and 
S J ) and the right ascension (a*) of the radio source cannot be separated in a common 
adjustment, since they are appeared in the partials through (Qi - a*) and (Q 2 -cx*). This 
fact indicates that a space-to-space VLB I time delay observable is insensitive in the location 
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of the true-of-date coordinate system where the coordinate difference components 
(coordinates of the endpoints) of a space-to-space baseline are given. Because of this 
linear interrelationship between these parameters, a new estimable parameter, namely Aa 0 ^ 
given by Eq. (41) must be introduced as unknown similarly to that of the cases of ground- 
to-ground and ground-to-space baseline VLBI observables. However, it can be shown that 
the equatorial components of the space-to-space baseline will be contaminated by the error 
doto of the initial reference right ascension in the following manner 

dAX^ 

dAY^ 

where dAx£ J and dAX^ are the differential changes in the equatorial components of the 
space-to-space baseline B IJ . 



1 

o- 

t> 
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0 

doto ’ 
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< 
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. dAY? . 


— dcto 

0 
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> 

l 


(179) 


3.6.2 Time Delay Rate Model. Differentiating Eq. (157) with respect to 


time we get 



Xi-X{ r 

y[-y{ 

zlc-z { L 


cosS^ COSOC/ 

cosS^ sina^ 
sin8/ 


+ cAC^ 


(180) 


where x£, Yk, andXj[, Yj[, Zj[ are the velocity components of the orbiting 
radiotelescope S 1 and S J , respectively, in the true-of-date inertial reference frame at the 
epoch tk . The model Eq. (180) contains nine unknown parameters 

{xL yL Z i xl y{, Z l a/, 5/, AC^ ) (181) 

to be determined. Taking the differential of dSr with respect to the parameters, 
d& = X F p dP p , 

P (182) 

where the F p ’s are the required partial derivatives of the time delay with respect to the 
parameters of interest, indexed by p . They are as follows: 


(183) 


Fxk = - cos 8; cos at 
Fyj^ = - cos 8/ sin a* 
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(184) 



(185) 


F^ k = - sin 8/ 

= cos 6/ cos a* 
Fy^ = cos 5/ sin a/ 
F^ = sin 5/ 


F a/ = cos 5/ [(Xk - Xk) sin a, - (Yk - yJ) cos a. 


F 5/ = sin 5/ (Xk - Xk) cos a t + (Yk - Yk) sin aj - (Z* - Z*) cos 5/ 


Fac? = c 

Irq 


(186) 

(187) 

(188) 

(189) 

(190) 

(191) 


Eqs. (183) - (188) indicate a complete linear dependence 

(Fj^i = -F^J , Fyi = -FyJ , F^ = -F^J) between the corresponding velocity components 
of the orbiting radiotelescope and S^, respectively. Since they cannot be separated at 
all, the velocity component differences of the orbiting radiotelescope S 1 and S J 

AXk = Xk — Xk (192) 

AYk = Yk - Yk (193) 

AZ?=z£-z£ (194) 


are, therefore, introduced as new combined unknowns with the partial derivatives 


F AX k U=F 4 


F AZ|? = F ZL 


(195) 

(196) 

(197) 


Instead of satellite’s Cartesian coordinates, Keplerian orbital parameters should be 
used in practice. For this purpose, the partials given in Eqs. (132) - (137) are useful. 
Using this parametrization, it can be shown that the ascending nodes (Oj and Q 2 ) of tho 
orbiting radiotelescopes and and the right ascension (<x g) of radio source are linearly 
interrelated, that is, they cannot be separated in a common adjustment. This inseparable 
situation for the parameters Qj, Q 2 and a £ indicates that the space-to- space baseline time 
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delay rate observable is also insensitive to the orientation of the true-of-date inertial 
reference system in right ascension. 

In order to circumvent to estimability problem, the right ascension of one 
(reference) source may be constrained to an adopted value and will not be included in the 
parameter set. Therefore, the corresponding estimable parameters are the right ascension 
differences Aa 0 * similarly to that of the previous observation models. However, the 
equatorial velocity components of the space-to- space baseline will be contaminated by the 
error d(Xo of the initial reference right ascension in the following manner: 

dAXk 1 _ [ dAX? I r 0 ddo 1 U*k ] , 19g , 

dAY'“ dAYj? J [-dOo 0 AY” 

where d AX^ and d AYj^ are the differential changes in the equatorial velocity components 
of the space-to-space baseline B IJ . 

3.7 Differential VLBI Observations 

Besides the usual VLBI observables (time delay, delay rate and phase delay), 
another very important tracking method is the differential VLBI (A VLB I) which can be 
applied to the determination of angular separation and its time derivative between a 
transponder on the orbiting radio observatory and natural extragalactic radio sources. As 
the name suggests, AVLBI is a measurement of the differential position between the 
satellite and a radio source; the measurement is made using VLBI and treating the satellite 
as a radio source. For an orbiting radiotelescope having a transponder on board, AVLBI 
observations from the ground of both the satellite and extragalactic radio sources in 
neighboring parts of the sky can serve to orient the satellite orbit with respect to the inertial 
frame formed by these sources. Since the positions of the natural radio sources are known 
with an accuracy of 0."001 as a result of astrometric programs, therefore the relative 
satellite coordinates are measurable with comparable accuracy if the satellite passes close to, 
and slowly by, one of the natural sources that constitute the inertial frame. 

In the AVLBI technique signals received simultaneously from a satellite and a radio 
source located close together in the sky are compared. The difference between the phases 
of the signals received simultaneously from a radio source at the two ground stations, 
usually called the interferometric “fringe phase”, is 
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(tk) — B ij Ok) • 6/ » 


(199) 


where (0 is the angular frequency of the transmitted signal, c is the speed of light, and 
By • e* is the pathlength difference [King, 1975]. 

If an orbiting transmitter and a natural radio source are observed simultaneously, 
the differential fringe phase is given by 

(tk) = <t>ij* (tk) - <t>ij Ok) = ^ B ij (tk) • (e* - ej;) (200) 

where we have supposed that the two transmitter frequencies are the same. Because 
AVLBI involves taking differences not only between receiving points but also between 
transmitting points, it follows that any potential source of error will cancel if it is common 
either to all receivers or to all transmitters. 

Expressing the small term (e*-e£) in terms of the small differential right 
ascension, Aa^ and differential declination, A 5^ between the directions to a natural radio 
source and an orbiting transmitter, we can write 

al=a/ + Aab (201a) 

§k = 8* + A8k* (201b) 

Whenever the orbiting radiotelescope passes close to the direction of an 
extragalactic radio source, AVLBI may be useful to determine the ground station — 
orbiting radiotelescope direction relative to that of the source to about 0!'001. Such 
measurements could be used to determine precisely the orientation of the radio observatory 
orbits with respect to an inertial frame and to improve the estimates of the Earth’s gravity 
model parameters. 
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4 . ESTIMATION PROCESS AND SOME NUMERICAL TESTS 


Some numerical computations were performed partly to check the mathematical 
models derived for space VLBI observables and partly to explore the overall applicability of 
space VLBI. In the following, first we discuss the datum definition problem in a general 
sense. This is followed by a short description of the numerical adjustment. Finally, results 
of the test computations are discussed. 


4.1 Datum Definition for the Space/Ground VLBI Network 

The general theory of network design including the geodetic datum problem and 
transformations is extensively discussed in the contributions, e.g., in Grafarend and Sanso 
(1985) and Teunissen (1985). In the following, the theoretical basis of the datum problem 
for a combined space/ground geodetic VLBI network will be shortly described. 

The parameters in Eq. (82) we are looking for cannot be derived uniquely from the 
space VLBI time delay measurements. This typical situation in a combined space/ground 
VLBI geodetic network can be characterized by a rank defect of the design matrix C within 
the Gauss-Markov Model, see [Schaffrin, 1985 and 1990], 

y - e = E(y) = Cx; o(C) = n x m , rk(C) = qcmcn , (202a) 

D(y) = D(e) = P -1 positive-definite, (202b) 

leading to the singular normal equations 


Nx = b ; N := C T PC , rk N = q<m , b := C T Py 


(203) 


Here y is the n x 1 observation vector, e the n x 1 error vector, C the n x m coefficient 
matrix (positive-semidefinite), x the m x 1 parameter vector, L ec the n x n positive- 
definite variance-covariance matrix of the observations, and P the respective n x n weight 
matrix (as well positive-definite). E denotes “expectation” and D “dispersion”. 

Eq. (203) yields a whole class of different estimates x 
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X = (C T PC)- C T Py = (C T PC) r - s C T Py 


(204) 


with accordingly different variance-covariance matrices 

D(x) = Nr s positive-semidefinite , (205) 

out of the class of reflexive symmetric g-inverses of N . Note that 

Cx = C(C T PC)-C T Py (206) 

is Best Linear Uniformly Unbiased Estimation (BLUUE) being unique for any choice of 
the g-inverse (C T PC)~ . 

From the rank deficiency of matrix C in (202a) [or (83)] follows that the unknowns x 
cannot be determined uniquely, even if y e R(C) . Thus, the information contained in y 
is not sufficient to determine x uniquely. Note that the question of how to define vector x 
so that it is estimable from y involves the so-called datum-problem. The datum-problem 
deals with coordinate system definition, since coordinates are not estimable quantities from 
observations alone. 

One can overcome this problem by adding the minimum information needed to 
determine x uniquely. Therefore, it is sufficient to introduce r additional minimal 
constraints 

Kx = Ko ; o(K) = r X m , rk(K) = r := m-q (207) 

such that the space conditions 

R(K T ) u R(C T ) = IR m (208a) 

and 

R(K T ) n R(C t ) = {01 (208b) 

hold, where R( ) denotes the range space of the corresponding matrices. From Eqs. (208a) 
and (208b) it follows that 


R(K t ) © R(C t ) = IR 111 


(209a) 



or 


R(K t ) © L(C) = R m , (209b) 

where L(C) is the nullspace of the design matrix C. (Note that the concepts of nullspace 
and range space are very important for singular matrices.) 

The famous restriction 

r = m-q (210) 

is for the number of minimal constraints which are needed to overcome the rank deficiency 
of the matrix C. Eq. (208b) says that the rows of the matrices C and K ought to be 
complementary (but not necessarily orthogonal). An equivalent condition with Eq. (208a) 
the 


N := N + K t K regular , 


( 211 ) 


is fulfilled. In this case, we obtain the estimated parameters x and their dispersion matrix 
D( x ) from the extended normal equations 


N 

K t ' 


X 


b ' 

K 

0 


_X_ 


.Ko. 


( 212 ) 


where X denotes some r x 1 vector of “Langrange multipliers.” 
In practice, special datum constraints 
Gx = 0; o(G) = r x m , rk(G) = r := m-q 

can be obtained through the orthogonality relation 
CG t = 0 . 


(213) 


(214) 


The rows of matrix G yield a basis of the nullspace L(C) so that condition Eq. (209a) then 
reads as 


R(G T )i R(C t ) = R m 


(215) 
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i.e., the range spaces R(C^) and R(G^) are perpendicular to each other, thereby 
automatically implying complementarity. 


If 


K := G 


(216) 


is available for an arbitrary, but fixed r x m matrix G with rk(G) = r and CG T - 0, then the 
“free net adjustment” leads to the estimates 

x = N + b (217a) 

D(x) = Oe 2 N + , (217b) 


where N + denotes the pseudo-inverse of N — C^PC . N + is found directly by 

N+ = (N + G t G) _1 - G t (GG t GG t )- 1 G , (218) 

which is independent of the choice of the matrix G. 

If G matrix is taken, the so-called minimum norm solution is found, which has a 
minimal value for the trace of the variance-covariance matrix. 

Adding the minimal constraints Eqs. (207) or (213) to the model Eq. (202a) removes 
the datum rank deficiency of the model. Applied to our combined space/ground VLBI 
geodetic network it means, that by adding the either minimal constraints the datum is 
established. 

Note that if we collect the eigenvectors belonging to the r-multiple eigenvalue 0 of 
the normal equation matrix C9C, we may even obtain an orthogonal basis (but it is not a 
necessary requirement). Therefore, we may apply a formal Singular Value Decomposition 
(SVD) to the positive semi-definite matrix C^C. Since C^C is a real symmetric semi- 
positive definite matrix, it follows that its eigenvalues are non-negative. Denoting these 
eigenvalues by of , i = 1,2,3,..., m, we can arrange them in decreasing order that 

Cl > a 2 > ... ^ <7 q > 0 = dq+1 = ... = (219) 
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For C T C, there exist q different positive eigenvalues and a r-fold (r = m-g) 
eigenvalue zero. Associated with these eigenvalues exist m different mutual orthogonal 
eigenvectors. 

The corresponding orthogonal eigenvectors are denoted by (vj, V 2 , ..., v m ) and we 
separate them into 

V 0 = (vi, V 2 , ..., v r ) (220a) 

Vi = (v q+ i, v q+ 2 , ..., v m ). (220b) 

With Aq 72 = diag (d, a 2 , ..., a q ) we thus have 

CTCVo^VoAq ( 221) 

and 

C T C Vi =0 (222a) 

or 

CVi = 0 (222b) 

Since the coordinate system defining subspace R(K T ) only needs to fulfill the 
complementarity condition Eq. (209b), it follows that there are many more coordinate 
frames possible, giving the same solutions in principle. Therefore, a relation exists 
between these solutions which is called S-transformation. For instance, a very special 
datum choice is given with Eq. (216). The condition Eq. (209b) is for practical purposes 
most manageable, since 

R(Vi)=L(C) (223) 

is given by the linearized S-transformation [Teunissen, 1985] and K (or K := G) is 
chosen in order to define the coordinate system. 

4.2 Numerical Adjustment 

In the mathematical models of space VLBI time delay and delay rate observables 
different type of parameters are included. However, in practice we are primarily interested 
in some of them, that is, e.g., earth rotation parameters (£, T|, k) are of most interest for 
geodynamical studies, or in order to improve the Earth’s gravity field, only the corrections 
to the Keplerian orbital elements are to be determined. Furthermore, supposedly we have 
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stochastic prior information on certain parameters. Therefore, a solution for a group of 
parameters in primary interest is often sought from a partitioned system of normal 
equations of the following form [Schaffrin, 1990] 


Nn Ni 2 

Xl 


CT 

N 2 1 N 22 . 

. x 2. 


1 

n 

NH 


(224) 


where 

C = [Cl C 2 ] (225) 

is the partitioning of the design matrix C of the space VLBI network and 

x = [xi x 2 ] (226) 

is the partitioning of the unknown parameters into a sub-vector xj with stochastic prior 
information and a sub-vector x 2 of the parameters of interest without prior information. 


The inversion of normal equation system (224) is 


Nn 

N 2 i 


N12 

n 22 J 


i-i 


sr 1 

-n 22 n 2 i sr 1 


-sr 1 Ni 2 r4 

N 2 i N 21 sr 1 Ni 2 N^ + N^ 


(227) 


where 

Si = Nh-Ni 2 N^N 21 (228) 

is the so-called “1st Schur complement”. Consequently, the estimated dispersion of the 
estimated parameters x 2 of interest is given by 

D {x 2 } = Oq (NjJ + N~ 2 \ N 21 Sr 1 N !2 N^) • (229) 

However, the system of normal equations cannot be inverted because of the datum rank 
defect. In practice the datum rank defect may be overcome by using certain conventions for 
the coordinate system definition, that is, e.g., by the use of a minimum norm solution (see 
the previous paragraph). In that case, in order to get the dispersion of the estimated 
parameters x 2 of interest, the pseudoinverse N + of N = C T PC given by Eq. (218) has to be 
used in partitioning in the practical computations. 
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4.3 Results of Test Computations 

In Section 3.5, the analytical partials and linear dependency are derived for both space 
VLB I observables. In the case of time delay observable, the rank defect of the normal 
equation matrix N = C T C is discussed in detail. Numerical checking of rank defect by 
evaluating numerical values of partials has been done in the frame of test computations. 
For this purpose we computed the normal equation matrix N = C^C . Eigenvalue analysis 
was used to check the rank defect of N = C T C . 

In order to derive the normal equation system matrix N=C T C numerically, generated 
data was used within a simulated ground/space VLBI network. For this purpose, 
minimum numbers were derived for participating stations, observed satellite positions and 
radio sources. The basic configuration consists of one satellite orbit for VSOP mission, 
three ground-based VLBI stations with three different radio source positions. 

The stations used in test computations are: Crimea, Jodrel 2 and Ovro 130. Their 
coordinates are given in Table 2. The used source positions as published in [Carter et al., 
1989] are collected in Table 3. 

In order to compute satellite spatial coordinates, initial orbital parameters of the VSOP 
satellite were used. A standard orbit simulation parameter set for VSOP is given in Table 
4. On the basis of these parameters, simulated orbit was generated for date of 01/01/1996 
using a space VLBI related program package developed at the Satellite Geodetic 
Observatory in Hungary and made available by I Fejes. After computation of satellite 
Cartesian position and velocity component, a set of “approximate” Keplerian orbital 
parameters was generated to use it for the numerical evaluation of analytical partials. 

The distribution of the three ground stations and the subsatellite track of the VSOP is 
shown in Fig. 17. An overall view of the stations and VSOP configuration is represented 
in Fig. 18. 

Simulated earth rotation parameters (^, T|, k) for epoch MJD50083 (01/01/1996) 
were calculated by using the prediction formulas given in [IERS Bulletin-A, 1990]. Their 
numerical values of ^ = -Q" 142486, T) = Q" 193437 and K = -0.17271 ms were used. 
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Satellite: VSOP 

Epoch: start 960101 00.00 

Stop 960102 00.00 


Fig. 17 Subsatellite track line of VSOP and Location of the 
Ground-Based VLBI Stations Used in the Test Computations 


Source: 1803+784 

ST: 56792d 2 h 24 m 

UT: 50083d 19 h 41 m 



SATELLITE : VSOP 

Fig. 18 Stations and VSOP Configuration Overall View 
from the Radio Source Direction 
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Table 2 

Station Coordinates Used in Test Computations 


Station 

X (m) 

Y (m) 

Z(m) 

Crimea 

378 5227.20 

255 1211.80 

443 9806.93 

Jodrel 2 

382 2842.66 

-15 3800.13 

508 6287.22 

Ovro 130 

-240 9626.30 

-447 8405.30 

383 8606.70 


Table 3 

Source Positions Used in Test Computations 

Source 

Right Ascension 

Declination of Arc 


H M S 

O i ft 

0212 + 735 

2 17 30.813210 

73 49 32.62230 

1641 + 399 

16 42 58.809930 

39 48 36.99380 

1803 + 784 

18 0 45.683850 

78 28 4.01790 


Table 4 

Initial Orbit Parameter Set for VSOP Satellite 


Keplerian Orbital Elements of VSOP Satellite 

Epoch of orbital elements 

47894.00000 [MJD] 

Semi major axis 

16878.00000 km 

Eccentricity 

0.56300 

Inclination 

31.000 deg 

Argument of perigee 

0.0000 deg 

Longitude of ascending node 

90.0000 deg 

Mean anomaly 

0.0000 deg 


Table 5 

Observations and Unknown Parameters of the Basic 
Configuration for Test Computations 


Observations 

Estimated Parameters 

12x3 time delavs 
36 

3x3 station coordinates 
3x2 source coordinates 
1x6 Keplerian elements 
1x3 ERP parameters 
3x2 clock parameters 
30 
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The Greenwich Apparent Sidereal Time (GAST = ©k) was computed from observed 
Universal Time UT1 data in MJD for the observation time epochs t k using formula given in 
[Obenson, 1970; p.25]. 

The relevant numbers which show the basic configuration for the test computation are 
shown in Table 5. On the basis of this configuration, there are 30 unknown parameters. 
Since we have three delay observations at each epoch, we need observations at ten different 
epoch for a unique solution. Instead of this minimum requirement, 12 satellite points were 
included to form our basic configuration for a space/ground VLBI network. Four satellite 
positions belong to one radio source observation period. 

The input metric data (station coordinates and semi-major axis) were used in Mm 
unit, and the angle type data unit was in radian. The speed of light was used in unit of 
Mm/msec. 

Using this generated data set, the partials of the time delay and delay rate with respect 
to the parameters involved have been calculated. Then the design matrix C for the time 
delay observable was found, and the elements of the normal equation system matrix N = 
C T C was computed. The dimension of the design matrix C and the normal equation system 
matrix N = C T C is 36x30 and 30x30 respectively. 

Two cases were considered in the computations concerning the orbit length: 1) short- 
arc (1/5 part of a full orbit), and 2) full orbit was used. In the first case, satellite positions 
were taken every six minutes, and for the second case the integration time period was 30 
mi nutes. 


It is not surprising that the normal equations derived from observation equations 
results in a singular system. This is confirmed by an eigenvalue analysis and using LU and 
Cholesky decomposition. 

The normal matrix N = C T C derived from the simulated data possesses a peculiar 
eigenvalue spectrum. The spectrum resembles that of a poorly conditioned system of 
equations with corresponding eigenvalues ranging from 10~^ to 10 + ^ in the case of Short- 
arc” configuration. In the “full orbit” situation, the eigenvalue spectra is in the range of 
10-11 to 10+ 4 . As a result, the condition number AmaxAminX representing a measure for 
the numerical stability of the system, is of the order of 10^ and 10^ respectively. In an 
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ideal situation, the condition number should lie in the range from 1 to 10 8 . This is not our 
case; that is, our normal equation system is poorly conditioned in both geometric 
configuration of the used space/ground VLBI network. 

The only cause of the poor condition number can be the observation equations 
themselves and results in an unstable normal equation matrix having approximately five 
order of magnitude difference in the values of its diagonal elements. The eigenvalue 
spectrum of the normal equation matrix can be improved by transformation of the design 
matrix C. It can be done with multiplication of the design matrix by a diagonal matrix P 
resulting a normal matrix N = C C having diagonal elements approximately the same 
order of magnitude in values. A multiplication of the design matrix C by a diagonal matrix 
P gives 

C = CP (230) 

and the normal equation system matrix becomes 

N = P r C T CP. (231) 

The inverse of this normal matrix N is (with P 1 = P) 

N _1 = (P C T C P) = P _1 (CTC)- 1 P _1 . 


The normal matrix N = C C , found after the transformation of the design matrix C, 
shows an eigenvalue spectrum having better properties. The results of this analysis are 
presented in Table 6. The computed eigenvalues of the normal matrix N show nine small 
eigenvalues indicating a considerable jump between the fourth and fifth eigenvalue, and in 
addition a system which is still poorly conditioned. The first four small eigenvalues may 
be due to true singularity. The eigenvalue spectra of the “full-orbit” network configuration 
is improved. It clearly shows four small eigenvalues which may be considered as zero 
with respect to the other 26 eigenvalues. The latter eigenvalues are in the range of 3.52 to 
3.67 x 10 8 . The four zero eigenvalues clarify our results for linear dependencies derived 
analytically. 
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Table 6 

Eigenvalues of the Normal Matrices in Case 
“Short-arc” and “Full Orbit” 


| No. | “Short-arc” 

“Full Orbit” | 

1 

8.27 x 10- 8 

2.88 x 10-' 1 

2 

2.37 x 10- 7 1 

3.24 x KH 

3 

4.42 x 10- 7 1 

5.33 x 10~ ; 

4 

8.14 x 10- 7 1 

7.01 x 10" 7 

5 

1.84 x 10- 4 

3.52 

6 

1.59 x 10-3 

9.18 

7 

6.85 x 1(H 1 

6.24 x 10 

8 

2.91 x W 1 1 

1.89 x 10 1 

9 

9.67 x 10-^ 

7.47 x 10 2 

10 

3.99 

6.13 x W 

11 

1.98 x 10 

1.43 x 10 4 

12 

9.60 x 10 

3.36 x 10 4 




30 

4.58 x 10# ~ 

! 3.67 x 10# 


The eigen vectors belonging to the four zero eigenvalues are the nullspace base 
vectors of the design matrix. They can be used to form the column vectors of a G matrix to 
derive the pseudo-inverse N+ in Eq. (218) for a minimum norm solution. Therefore, 
another numerical test is to check the orthogonality of these eigen vectors with respect to 
the column vectors of the design matrix, that is, to check the numerical validity of Eqs. 
(222a) and (222b). The results obtained for the elements of matrix multiplication expressed 
by Eqs. (222a) and (222b) are of the order of 10 -6 to 10 -4 . 

As discussed in section 3.5.1 singularity of the normal equations is caused by linear 
dependency of the column vectors in the design matrix. In order to confirm our results 
given by Eqs. (108) and (109) in another way, the so-called LU-decomposition and 
Cholesky decomposition of the normal equation system matrix N = C C have been 
performed as an independent procedure. For LU-decomposition a routine in [Press et al., 
1986] has been used. The Cholesky reduction algorithm is given by Stewart (1973; p. 
142). The parameters to be solved for are involved in the normal equation system in the 
following order: 

{Xi, Yi, Zi, X 2 , Y 2 , Z 2 , X 3 , Y 3 , Z 3 , a, e, i, CD, Q, M, k, 

<xi, 8i, a 2 , 8 2 , a 3 , 8 3 , AC<j rl , ACi rl , ACo r2 , AC{ r2 , ACo r3 , ACi r3 ) (223) 
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Both types of decomposition of normal matrix show small values for the diagonal elements 
corresponding to the earth rotation parameters and the right ascension of the third radio 
sources. In the case of LU-decomposition of C C , the small values in question are in the 
range of lfK 7 to 10 -3 . The corresponding values from Cholesky decomposition are of the 
order of 10 -3 . They may be considered as zero with respect to the other diagonal elements 
in both decomposed forms of normal matrix. A zero value means complete dependence on 
(a linear combination of) other unknown parameters. 

A simple checking for linear dependencies has also been performed by calculation 
of the numerical values of the Eqs. (108) - (109) and Eqs. (149) - (150). The validity of 
Eqs. (Ill) -(113) and Eqs. (151) - (153) derived for critical configuration study has been 
checked numerically as well. 
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5. SUMMARY AND RECOMMENDATIONS 


Investigations presented in this study are summarized as follows: 

(1) Qualitative expressions have been developed to demonstrate the relationship 
between the space VLBI observables (time delay and delay rate) and the solve for 
parameters, and to explore estimability problems inherent in the space VLBI system. 
Instead of Cartesian satellite coordinates, Keplerian orbital elements were introduced into 
the observation equations and partial derivatives of space VLBI observables with respect to 
the parameters of interest were derived. A sensitivity analysis using these partials has been 
carried out. The mathematical models for ground-to-ground baseline VLBI observables are 
also given for comparison. Simplified models were used for the space-to-space VLBI 
observation equations. 

(2) Linear dependencies between partial derivatives have been derived analytically. 
Their expressions for both space VLBI (ground-to- space baseline time delay and delay rate) 
observables are formally similar. Four independent linear dependencies were found 
indicating that the number of the rank defect of the corresponding normal equation matrix is 
four. It has been shown through analytical derivation that the linear dependencies between 
partials of the ground-to-space baseline VLBI observables are formally very similar to ones 
of the ground-to-ground VLBI observables. However, in the case of both ground-to-space 
baseline VLBI observables the expression of linear dependency between earth rotation 
parameter k and the right ascension of radio source a e were extended by a new term for 
ascending node Q of the satellite orbit indicating that these three parameters are 
nonseparable in a common adjustment. 

(3) A geometric interpretation of datum rank defect (analytically derived linear 
dependencies) in the space VLBI time delay observation equations has been presented. The 
rank defects expressed by linear dependencies are confirmed using different numerical 
methods with simulated data set in the frame of test computations. The eigenvalue analysis 
and both decomposition (Cholesky and LU of a lower and an upper triangular matrix) of 
the normal equation matrix indicate that the number of datum rank defects is four 
concerning the mathematical model used for the space VLBI time delay observable. 
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the normal equation matrix indicate that the number of datum rank defects is four 
concerning the mathematical model used for the space VLBI time delay observable. 

(4) In searching linear dependencies, analytical expressions were derived to study 
critical configuration in the case of both ground-based and space VLBI. They are useful in 
observation design work. A check of these relations was performed by test computation. 

(5) Datum definition problem for space VLBI network is theoretically described 
within the framework of rank defect Gauss-Markov Model. Proper numerical adjustment 
procedure is given which can be used for covariance analysis of the parameters of geodetic 
interest. 

In order to use successfully and efficiently the space VLBI observations in geodesy 
and geodynamics, several additional studies need to be made in order to investigate in full 
detail the potential advantages and the possible problems of the space VLBI technique. 
There have been no simulation studies so far to indicate the expected accuracies of the space 
VLBI system in the determination of the parameters of astrometric and geodetic- 
geodynamic interest. Therefore, a feasibility study should be undertaken in order to show 
whether or not space VLBI can provide real improvements in the potential application in 
geodesy and geodynamics. Hence, in order to give a feeling of what one should expect 
from such a system, a detailed sensitivity analysis needs to be carried out in order to 
estimate the accuracy achievable in determining various parameters. More emphasis should 
be given to the determination of earth rotation parameters for reference frame studies and to 
the estimation of orbit parameters to improve the gravity field determination. On the basis 
of this work, we would recommend that: 

a) The mathematical model be augmented to include additional effects like variation 
of Keplerian orbital elements. Using the extended mathematical model, the rank defect 
analysis should be repeated. 

b) Simulations in order to determine optimum global station network and orbit 
configuration, and accuracy estimates for the various parameters are needed. These would 
require extended simulation studies and error analysis using data from a realistic 
observational pattern. 
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c) Another area which needs to be developed in the near future is the software 
development for different tasks involved in space VLBI geodetic data analysis. 

d) Such investigations should provide sufficient background for the inclusion of 
space VLBI observable in geodetic data processing programs (e.g., GEODYN). 
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APPENDIX A 


Keplerian — Cartesian Transformations 


First, the transformation between Keplerian orbital elements (a, e, i, co, £2, M} and 
Cartesian orbital (rectangular) coordinates {X, Y, Z, X, Y, Z) is considered. This is 
followed by the description of partial derivatives of the Cartesian orbital coordinates with 
respect to the Keplerian elements and partial derivatives of the Keplerian elements with 
respect to the Cartesian elements. 


A 1 . Kepler Transformation of Coordinate and Velocity Type 

Transformation of Kepler elements {a, e, i, to, £1, M) into inertial Cartesian 
position vector r : 



X 


cos Q. cos u - sin £2 sin u cos i 

r - 

Y 

= r 

sin Q cos u + cos Q sin u cos i 


_Z_ 


sin u sin i J 


(A.l) 


where 


r = a (1 -e cos E) 
u = to + f 

tg f/2 = V ( 1 +e)/(l-e) tg E/2 
M = E - e sin E 

sin f = V 1 - e 2 sin E/( 1 - e cos E) 
cos f = (cos E - e)/(l - e cos E) 

E is the eccentric anomaly. 


(A. 2) 
(A. 3) 
(A. 4) 
(A. 5) 
(A. 6) 
(A. 7) 
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Transformation of Kepler elements {a, e, i, co, £2, M} into inertial Cartesian velocity vector 

r: 



X 


cos £2 cos u - sin £2 sin u cos i 

r = 

Y 

= v r 

sin £2 cos u + cos £2 sin u cos i 


_Z_ 


sin u sin i 


- v 


u 


cos £2 sin u + sin £2 cos u cos i 
sin £2 sin u - cos £2 cos u cos i 


- cos u sin i 


where 

v r = Yp./p e sin f 
v u = Y|i/p (1 +ecosf) 
p = a (1 -e 2 ) 

I x is the gravitational constant times the mass of the Earth. 


(A. 8) 


(A.9) 

(A. 10) 
(A. 11) 


An equivalent conversion from Keplerian elements to inertial Cartesian position and 
velocity vectors (r and r): 

a cos E - a e 
a Y 1 - e 2 sin E 
0 

- sin E 

Y 1 -e 2 cos E 
0 

Rxq = R 3 (-Q) • RlH) ■ Rs (-«) 

r = RxqR 
r = Rxqq 

where n is the mean motion. 

A 2 . Conversion from Inertial Cartesian Position and Velocity Vectors to Keplerian 
Elements (Inverse Kepler Transformation) 

Compute radial distance r and tangential velocity v : 

r = | r| = Yx 2 + Y 2 + Z 2 (A. 17) 


q = ^ 

(1 -e cos E) 


q = 


a (1 -e 2 ) 
(1 + e cos f) 


r ”1 


cos f 


sin f 

= 

0 . 
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V = | r| = tx 2 + Y 2 + Z 2 


(A. 18) 


Compute angular momentum vector h : 

{= 1 } 

h = lfh5+ hT+hJ 

Compute inclination angle i : 
i = cos 1 (V| 


h = r x r = 


YZ-ZY 

ZX-XZ 

XY-YX 


Compute ellipse semi-major axis a 
v 2 = GM 


2_1 

r a 


a = 


2 v 2 
r GM 


Compute ellipse eccentricity e : 
h 2 = GM a (1 -e 2 ) 

lV2 


e = 


h 2 

aGM 


Compute longitude of ascending node Q : 

n=tanI &) 

Compute true anomaly f : 

f = tan -1 f sin 0 , 

^cos f ) 

where 


(A. 19) 
(A. 20) 

(A.21) 

(A. 22) 
(A.23) 

(A. 24) 
(A. 25) 

(A. 26) 

(A. 27) 
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(A. 28) 


cos f = 


sin f = 


a(l-e 2 ) i 
er e 

a (1 -e 2 ) (r • f) 
reh 


Compute eccentric anomaly E : 

E = tan' 1 f sin ^ 1 
ycos E J 

where 

cos E = -G QS . f + e 
1 + e cos f 

. „ V 1 -e 2 sin f 

sm E = — — 

1 + e cos f 

Use Kepler’s equation 
M = E - e sin E 

to compute the mean anomaly M . 


Compute geocentric angle u : 

u = tan' 1 f sin u 
t^cos u 

where 

cos u _ X cos Q + Y sin Q 
r 

(Y cos Q - X sin Q) cos i + Z sin i 

sin u = 

r 

Compute argument of perigee co : 
to = u — f 


(A. 29) 

(A. 30) 

(A. 31) 
(A. 32) 

(A. 33) 


(A.34) 

(A. 35) 
(A.36) 

(A. 37) 
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A 3 . Partial Derivatives of Inertial Cartesian Coordinates and Velocity Components 
with respect to Keplerian Elements (Linear Kepler transformation) 

Elements of the linear Kepler transformation matrix of position type: 


ax 

da 


= (1 - e cos E) (cos Q cos u - sin Q sin u cos i) 


0Y 

— = (1 - e cos E) (sin Q cos u + cos £2 sin u cos i) 
da 


dZ 

— — = (1 - e cos E) sin u sin i 
da 


ax 

de 


= [a(e -cos E)/(l -e cos E)] (cos £2 cos u- sin Cl sin u cos i) 


- a sin f 


f , 1 - e cos E 

1 + r— 

1 -e 2 


(cos £2 sin u + sin £2 cos u cos i) 


— = [a(e - cos E)/(l - e cos E)] (cos £2 cos u + sin £2 sin u cos i) 
de 


a sin f 


( , 1 - e cos E^ 


1 + 


1 -e 2 


(sin £2 sin u - cos £2 cos u cos i) 


dZ 


— — = [a(e - cos E)/(l - e cos E)1 sin u sin i + a sin f 
de 


. „ f , 1 - e cos E ^ 


1 + 


1 -e 2 


cos u sin i 


ax 

dco 

dY 

do) 

dZ 

dco 

dX 

di 

dY 

di 


= - r (cos Q sin u + sin Q cos u cos i) 


= - r (sin Q sin u - cos Q cos u cos i) 


= r cos u sin i 


= r sin Q sin u sin i 


= - r cos Q sin u sin i 


(A. 38) 
(A. 39) 
(A. 40) 

(A. 41) 

(A. 42) 

(A.43) 
(A. 44) 
(A. 45) 
(A. 46) 
(A. 47) 
(A. 48) 
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(A.49) 


az 

ai 


= r sin u cos 1 


ax 

3£2 


= - r (sin £2 cos u + cos £2 sin u cos i) 


3Y 

= r (cos £2 cos u - sin Q sin u cos i) 

3£2 


az 

3£2 


= 0 


ax 

3M 


= [a e sin E/( 1 - e cos E)] (cos £2 cos u - sin £2 sin u cos i) - 


— [a sin f/sin E] (cos £2 sin u + sin £2 cos u cos i) 


= [a e sin E/(l - e cos E)] (sin £2 cos u + cos £2 sin u cos i) 

3M 

- [a sin f/sin E] (sin £2 sin u - cos £2 cos u cos i) 


^2 

- — = [a e sin E/(l - e cos E)] sin u sin i + [a sin f/sin E] cos u sin i 
3M 


Elements of the linear Kepler transformation matrix of velocity type: 


ax 


= - - 1 - (cos £2 cos u — sin £2 sin u cos i) + ^ (cos £2 sin u + sin £2 cos u cos i) 


3a 2a 


2a 


3Y 

3a 


— — = - (sin £2 cos u + cos £2 sin u cos i) + ^ (sin £2 sin u - cos £2 cos u cos i) 


2a 


2a 


az 

3a 


VjL - sin u sin i + ^ cos u sin i 


2a 


2a 


(A. 50) 
(A. 51) 
(A.52) 

(A. 53) 

(A. 54) 
(A. 55) 

(A. 56) 
(A.57) 
(A.58) 
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Bx 

Be 


- v, 1 (cos Q cos u - sin Q cos i sin u) 

e (1 -e cos E) 


BY 

Be 


BZ 

Be 


BX 

Bco 

BY 

Bco 

BZ 

Bco 


- v„ e + COS f ( cos Q sin u + sin Q cos i cos u) 

(1 -e 2 )(l + e cos f) 


= — v 


r e (1 -e cos E) 


- — (sin Q cos u + cos C2 sin u cos i) - 


v u e + cos f (sin q sin u — cos Q cos u cos i) 

(1 -e 2 )(l + e cos f) 


= — v 


1 


sin u sin l 


+ v u 


r e (1 -e cos E) 

e ± CQ S . f cos u s j n i 


(1 -e 2 )(l + e cos f) 


= - v r (cos Q sin u + sin Q. cos u cos i) - v u (cos Q cos u - sin Q sin u cos i) 


= - v r (sin Q sin u - cos Q cos u cos i) - v u (sin £2 cos u + cos f2 sin u cos i) 


= v r cos u sin l - v u sin u sin 


BX 

Bi 

BY 

Bi 

BZ 

Bi 


= v r sin Q sin u sin i + v u sin Q cos u sin i 


= - v r cos u sin u sin i - v u cos Q cos u sin i 


= v r sin u cos i + v u cos u cos i 


(A. 59) 


(A. 60) 

(A. 61) 

(A. 62) 
(A. 63) 
(A. 64) 

(A. 65) 
(A. 66) 
(A. 67) 
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= - v r (sin Q cos u + cos Q. sin u cos i) + v u (sin Q sin u - cos Q cos u cos i) 

3Q (A. 68) 

= v r (cos Q cos u - sin £2 sin u cos i) - v u (cos Q sin u + sin Q cos u cos i) 

(A. 69) 


az 

an 


= 0 


(A.70) 


= v r cos f — v u sin f q cos u - sin n sin u cos i) 
aM sin E (1 - e cos E) 


(A.71) 


3Y 

aM 


v r cos f - v u sin f 
sin E (1 — e cos E) 


(sin n cos u + cos n sin u cos i) 


(A.72) 


az 

aM 


v r cos f - v u sin f 
sin E (1 - e cos E) 


sin u sin i 


(A. 73) 


A 4 . Partial Derivatives of Keplerian Elements with respect to Inertial Cartesian 

Coordinates and Velocity Components (Inverse Linear Kepler Transformation) 

Notation: 

s = X, Y, Z, X, Y, Z respectively (A.74) 


Partial derivative of the radial component: 

5* 0, o, 0 respectively 

as r r r 


Partial derivative of the velocity component: 
= 0, 0, 0, ^ | respectively 


Partial derivative of the angular momentum h 

ah i f ah x , ah Y , ah z 

-r- = r- hx ^r— + hy ^r— + h Z — 
as h d s os as 


(A. 75) 


(A. 76) 


(A.77) 
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where 


~ = 0, Z, -Y, 0, -Z, Y 
ds 

respectively , 

(A. 78) 

^ = -Z, 0, X, Z, 0,-X 
ds 

respectively , 

(A. 79) 

^ = Y,-X, 0, -Y, X,0 
ds 

respectively . 

(A. 80) 


Partial derivative of (r • r) : 

r • r = XX + YY + ZZ (A 81) 

8 (l ~ ' r) = X, Y, Z, X, Y, Z respectively 

ds (A.82) 


Partial derivative of the inclination i : 


*=± 
3s h 2 


'hxh z ^ + h Y h z ^-(hi + h^l 

ds ds ds 

(h^+h ^) 1/2 


(A. 83) 


Partial derivative of the semi-major axis a : 


3a 

3s 


r 2 . 3r + 2 w_ > 

r 2 3s GM 3s 

J 


3a 

3s 


2a 2 
r 2 3s 


2a 2 

GM 



s = X, Y, Z 
s = X, Y, Z 


Partial derivative of the eccentricity e: 

3e _ 1 - e 2 3a 1 - e 2 3h 
3s 2ae 3 s he 3s 


(A. 84) 


(A. 85) 


(A. 86) 
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Partial derivative of the ascending node Q , : 


(A. 87) 


(A. 88) 


(A. 89) 


(A. 90) 
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Partial derivative of u : 


— = 1 -Tf- Y sin Q cos i cos u - X cos Q cos i cos u 
3s r 1L 

aa 


+ X sin Q sin u - Y cos Q sin u] 


a s 


3i 

+ [X sin Q sin i cos u - Y cos Q sin i cos u + Z cos i cos u] — 
L ds 

,ax 


+ [- sin Q cos i cos u - cos £2 sin u] 
+ [cos Q cos i cos u - sin Q sin u] 

r . 1 aZ 1 

+ [sin l cos uj — ] 

Partial derivative of the argument of perigee co : 

aco _ au af 

as 3s 3s 


as 

3Y 

3s 


(A.91) 


(A. 92) 


A5. Let the matrix K c representing the linear Kepler transformation matrix whose 
elements are the partial derivative of the Cartesian orbital elements with respect to the 
Keplerian elements: 


ax 

ax 

ax ax ax ax 

3a 

3e 

di d(D dQ dM 

3Y 

3Y 

3Y 3Y 3Y 3Y 

3a 

3e 

di 3 co 3Q 3M 

az 

az 

3Z 3Z 3Z 3Z 

3a 

3e 

3i 3 co 3£1 3M 

ax 

ax 

ax ax ax ax 

3a 

3e 

3i 3 co 3Q 3M 

3Y 

3Y 

3Y 3Y 3Y 3Y 

3a 

3e 

3i 3 co a O 3M 

az 

az 

az az az az 

3a 

3e 

3i 3 co 3Q 3M 


(A. 93) 
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Let the partial derivatives of the Keplerian orbital elements with respect to the 
Cartesian elements be the elements of the inverse linear Kepler transformation matrix K c in 
the following form: 


8a 0a 0a 0a 0a 0a^ 
0X 0Y 0Z ax 0Y 0Z 
0e 0e 0e 0e 0e 0e 

0x0Y0za^aV0^ 
0i 0i 0i 0i 0i 0i 
0x 0 y 0z ax at 0z 
0(0 00 ) 00 ) 00 ) 00 ) 00 ) 
0x 0Y0z a^at 0z 
an 0Q 0Q 0Q bq 0a 
0X 0Y 0Z ax 0Y 0Z 
aM 0M 0M 0M 0M 0M 

ax dr 0z ax by az 


(A.94) 


Then 


K C K C = K C K C = I . 


(A.95) 


where K c and K c are the Jacobians of transformation from Keplerian orbital elements into 
Cartesian coordinates respectively from Cartesian coordinates into Keplerian orbital 
elements. 

In order to be safe from the unavoidable typing errors in the relations given in this 
Appendix, the program listing of these relations is given in Appendix B for a comparison 
and practical use. The program was tested for numerical data. The output values of a test 
run are also given in Appendix B. 
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APPENDIX B 


Fortran Program DOKEPLER for Keplerian-Cartesian Transformations 


The program is written in Fortran 77 and consists, beside the main program unit, of 
two subroutines. The main program unit reads the input data, organizes the work of the 
subroutines and creates two new data files. The subroutines are as follows: 

1) Subroutine LEKEPTRA calculates the Cartesian position and velocity components 
of the satellite and their partial derivatives with respect to the Keplerian orbital elements, 
provided the Keplerian orbital elements (a, e, i, co, £2, M, and f,E) and the Earth’s 
gravitational constant, GM are known. The output values are the Cartesian satellite 
coordinates X,Y,Z and velocity components X,Y,Z as well as the elements of the 6 x 6 
Jacobian Kc of transformation from Keplerian orbital elements into Cartesian coordinates. 

2) Subroutine INKEPTRA calculates the Keplerian orbital elements (a, e, i, co, Q, 
M, and f,E) together with the elements of the 6x6 Jacobian K c of transformation from 
Cartesian coordinates into Keplerian orbital elements. The input data, which should be 
defined before calling the subroutine, are the rectangular components of position and 
velocity of the satellite and the Earth’s gravitational constant. 

The program calculates the elements of both multiplication matrices of the Jacobians, 
Kc and K c in order to check Eq. (A.95) in Appendix A. 

After the listing of the program and the subroutines, this Appendix contains numerical 
examples for Keplerian-Cartesian Transformations and for checking the multiplication 
products of the Jacobians. Instantaneous VSOP satellite position and velocity components 
respectively Keplerian orbital elements at epoch MJD 50083.000 (or GAST = 
105854.567125788) have been used. 
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* 

* 


* 


***** PROGRAM dokepler ***** * 

TO COMPUTE THE KEPLERIAN-CARTESIAN TRANSFORMATIONS AND THE * 
ELEMENTS OF THE CORRESPONDING LINEAR KEPLER TRANSFORMATION * 

MATRICES ( KC , KCI ) USING SUBROUTINES LIKEPTRA AND INKEPTRA * 

FINAL PROGRAM TESTED FOR DATA * 

★ ★★★★★★★★★*******‘**********‘***** , * r ** , ** r ********** r *************** * 


IMPLICIT REAL* 8 ( A-H, O-Z ) 

DOUBLE PRECISION KC ( 6 , 6 ) , KCI (6,6) 

DIMENSION UN (6,6), UNI (6,6) 

OPEN ( UNIT=1 , FILE= ' A : CDATACKEPLERCA . DAT ' ) 

OPEN(UNIT=2 , FILE= ' A : CDATACKEPLERCA . OUT ' ) 

OPEN ( UNIT=3 , FILE= ' A s CDATACKEPLERUN . OUT ' ) 

READ ( 1 , * ) AO , EO , XIO 
READ ( 1 , * ) OMO , OMEO , XMO 
READ ( 1 , * ) FO , EAO 
READ ( 1 , * ) XO , YO , ZO 
READ ( 1 , * ) XDO , YDO , ZDO 

CALL LIKEPTRA ( AO , EO , XIO , OMO , OMEO , XMO , FO , EAO , XS , YS , ZS , XSD , 
*YSD, ZSD,KC) 

WRITE ( 2 , 6 ) XS , YS , ZS 
WRITE (2,6) XSD , YSD , ZSD 

6 FORMAT ( 3 ( IX, F23 . 11 ) ) 

DO 1 1=1,6 

WRITE (2,7) (KC(I,J) ,J=1,3) 

WRITE (2,7) ( KC( I , J ) , J=4 , 6 ) 

1 CONTINUE 

7 FORMAT (3(1X,F23.13) ) 

CALL INKEPTRA ( XO , YO , Z 0 , XDO , YDO , ZDO , AS , ES , XIS , 

*OMS , OMES , XMS , FS , EAS , KCI ) 

WRITE { 2 , * ) AS , ES , XIS 
WRITE ( 2 , * ) OMS , OMES , XMS 
WRITE ( 2 , * ) FS , EAS 
DO 2 1=1,6 

WRITE (2,8) (KCI(I,Cn ,J=1,3) 

WRITE (2,8) ( KCI { I , J ) , J=4 , 6 ) 

2 CONTINUE 

8 FORMAT (3(1X,F23.17) ) 

CALL MATMULT ( KC , KCI , UN , 6 , 6 , 6 ) 

DO 3 1=1,6 

WRITE (3,4) ( UN ( I , J ) , J=1 , 3 ) 

WRITE (3,4) (UN(I, J) , J=4,6) 

3 CONTINUE 

4 FORMAT ( 3 ( IX, F23 . 20 ) ) 

CALL MATMULT ( KCI ,KC, UNI, 6 ,6, 6) 

DO 5 1=1,6 

WRITE (3,4) ( UNI ( I , J ) , J-l , 3 ) 

WRITE (3,4) (UNI(I,J),J-4,6) 

5 CONTINUE 
CLOSE ( 1 ) 

CLOSE ( 2 ) 

CLOSE ( 3 ) 

STOP 

END 

C$ INCLUDE AsCSVLBICLIKEPTRA 
C$ INCLUDE AsCSVLBICINKEPTRA 
C$ INCLUDE As CSVLBI<: MATMULT 




* ***** SUBROUTINE LIKEPTRA ***** * 

* TO- CONVERT THE KEPLERIAN ELEMENTS INTO SATELLITE'S CARTESIAN * 

* POSITION AND VELOCITY VECTORS AND TO COMPUTE THE ELEMENTS OF * 

* THE LINEAR KEPLER TRANSFORMATION MATRIX (KC) * 


* *************************************************************** * 
SUBROUTINE LIKEPTRA( A,E,XI,OM,OME,XM,F,EA,X,Y, Z,XD, YD, ZD,KC) 
IMPLICIT REAL *8 ( A-H, O-Z ) 

DOUBLE PRECISION KC(6,6) 

GM=398600 . 436D9 
C TO COMPUTE THE COEFFICIENTS 
U=OM+F 

COM=DCOS ( OME ) 

SOM=DSIN ( OME ) 

SU=DSIN(U) 

CU=DCOS ( U ) 

SI-DSIN(XI) 

CI=DCOS ( XI ) 

SF=DSIN ( F ) 

CF=DCOS ( F ) 

SE=DSIN(EA) 

CE=DCOS ( EA ) 

P=A* ( 1 . 0-E*E ) 

R=A* ( 1 . 0-E*CE ) 

VR=E*SF*DSQRT ( GM/P ) 

VU= ( 1 . 0+E*CF ) *DSQRT ( GM/P ) 

C TO COMPUTE SATELLITE COORDINATES 
X=R* ( COM*CU-SOM*SU*CI ) 

Y=R* ( SOM*CU+COM*SU*CI ) 

Z=R*SU*SI 

C TO COMPUTE SATELLITE VELOCITY COMPONENTS 

XD=VR* ( COM*CU-SOM*SU*CI ) -VU* ( COM*SU+SOM*CU*CI } 

YD=VR* ( SOM*CU+COM*SU*CI } -VU* ( SOM*SU-COM*CU*CI ) 
ZD=VR*SU*SI+VU*CU*SI 

C TO COMPUTE THE COEFFICIENTS OF THE PARTIALS 
Al= ( A* ( E-CE ) ) / ( 1 . 0-E*CE ) 

A2=l . 0+ ( 1 . 0-E*CE )/( 1 . 0-E*E ) 

A3= ( A*E*SE ) / ( 1 . 0-E*CE ) 

A4=-1.0/(E*(1.0-E*CE) ) 

A5= ( E+CF ) / ( (1.0-E*E)*(1.0+E*CF) ) 

A6= ( VR*CF-VU*SF ) / ( SE* ( 1 . 0-E*CE ) ) 

C TO COMPUTE THE PARTIALS of X 

KC ( 1 , 1 ) = ( 1 . 0-E*CE ) * ( COM*CU-SOM*SU*CI ) 

KC (1,2) =A1* ( COM*CU-SOM*SU*CI ) -A*SF*A2* (COM*SU+SOM*CU*CI ) 

KC ( 1 , 3 ) =R*SOM*SU*SI 

KC ( 1 , 4 ) =-R* ( COM*SU+SOM*CU*CI ) 

KC ( 1 , 5 ) =-R* ( SOM*CU+COM*SU*CI ) 

KC ( 1 , 6 ) =A3 * ( COM*CU-SOM*SU*CI ) - ( A*SF/SE ) * ( COM*SU+SOM*CU*CI ) 

C TO COMPUTE THE PARTIALS of Y 

KC ( 2 , 1 ) = ( 1 . 0-E*CE ) * ( SOM*CU+COM*SU*CI ) 

KC ( 2 , 2 ) =A1* ( SOM*CU+COM*SU*CI ) -A*SF*A2* ( SOM*SU-COM*CU*CI ) 

KC ( 2 , 3 ) =-R*COM*SU*SI 

KC ( 2 , 4 ) =-R* ( SOM*SU-COM*CU*CI ) 

KC (2,5) =R* { COM*CU-SOM*SU*CI ) 

KC ( 2 , 6 ) =A3 * ( SOM*CU+COM*SU*CI ) - ( A*SF/SE ) * ( SOM*SU-COM*CU*CI ) 

C TO COMPUTE THE PARTIALS of Z 
KC ( 3 , 1 ) = ( 1 . 0-E*CE ) *SU*SI 
KC (3,2) =A1*SU*SI+A*SF*A2*CU*SI 
KC (3,3) =R*SU*CI 
KC (3,4) =R*CU*SI 



KC( 3, 5 )=0 . 0 

KC ( 3 , 6 ) =A3*SU*SI+ ( A*SF/SE ) *CU*SI 
C TO COMPUTE THE PARTIALS of XD 

KC ( 4 , 1 ) = ( -VR* { COM*CU-SOM*SU*CI ) +VU* { COM*SU+SOM*CU*CI ) ) / ( 2 . 0*A ) 
KC ( 4 , 2 ) =VR*A4 * { COM*CU-SOM*CI *SU ) -VU*A5* ( COM*SU+SOM*CI *CU ) 

KC ( 4 , 3 ) =VR*SOM*SU*SI+VU*SOM*CU*SI 

KC ( 4 , 4 ) =-VR* ( COM*SU+SOM*CU*CI ) -VU* { COM*CU-SOM*SU*CI ) 

KC ( 4 , 5 ) =-VR* ( SOM*CU+COM*SU*CI ) +VU* ( SOM*SU-COM*CU*CI ) 

KC ( 4 , 6 ) =A6 * ( COM*CU-SOM*SU*CI ) 

C TO COMPUTE THE PARTIALS of YD 

KC ( 5 , 1 ) = { -VR* ( SOM*CU+COM*SU*CI ) +VU* ( SOM*SU-COM*CU*CI ) ) / ( 2 . 0*A ) 
KC ( 5 , 2 ) =VR*A4* ( SOM*CU+COM*CI*SU ) -VU*A5* ( SOM*SU-COM*CI*CU ) 

KC ( 5 , 3 ) =-VR*COM*SU*SI-VU*COM*CU*SI 

KC ( 5 , 4 ) =-VR* ( SOM*SU-COM*CU*CI ) -VU* ( SOM*CU+COM*SU*CI ) 

KC ( 5 , 5 ) =VR* ( COM*CU-SOM*SU*CI ) -VU* ( COM*SU+SOM*CU*CI ) 

KC ( 5 , 6 ) =A6 * ( SOM*CU+COM*SU*CI ) 

C TO COMPUTE THE PARTIALS of ZD 

KC ( 6 , 1 ) = ( -VR*SU*SI-VU*CU*SI ) / ( 2 . 0*A ) 

KC ( 6 , 2 ) =VR*A4*SU*SI+VU*A5*CU*SI 

KC ( 6 , 3 ) =VR*SU*CI+VU*CU*CI 

KC (6,4) =VR*CU*SI-VU*SU*SI 

KC( 6, 5 )=0 . 0 

KC (6,6) =A6*SU*SI 

RETURN 

END 


* **************************************************************** * 

* ***** SUBROUTINE INKEPTRA ***** * 

* TO CONVERT THE SATELLITE'S CARTESIAN POSITION AND VELOCITY * 

* VECTORS INTO KEPLERIAN ORBITAL ELEMENTS AND TO COMPUTE THE * 

* ELEMENTS OF THE INVERZ LINEAR KEPLER TRANSFORMATION MATRIX * 

* **************************************************************** * 

SUBROUTINE INKEPTRA ( X , Y , Z , XD , YD , ZD , A , E , XI , OM , OME , XM , F , EA , KCI ) 
IMPLICIT REAL* 8 ( A-H,0-Z ) 

DOUBLE PRECISION KCI (6, 6) 

GM=398600 . 436D9 
R=DSQRT(X**2+Y**2+Z**2) 

V=DSQRT ( XD* * 2+YD* * 2+ZD* * 2 ) 

RRD=X*XD+Y*YD+Z*ZD 
C MOMENTUM VECTOR 

HX=Y*ZD-Z*YD 

HY=Z*XD-X*ZD 

HZ=X*YD-Y*XD 

H=DSQRT{HX**2+HY**2+HZ**2 ) 

C INCLINATION (i) 

XI=DACOS ( HZ/H ) 

C SEMI -MAJOR AXIS (a) 

A= ( R*GM ) / ( 2 . DO *GM-R* V*V ) 

C ECCENTRICITY (e) 

E=DSQRT{ 1 . DO- ( H*H ) / ( A*GM) ) 

C LONGITUDE OF ASCENDING NODE (OMEGA) 

OME =D AT AN 2 ( HX, -HY ) 

C TRUE ANOMALY (f) 

COSF= ( A* ( 1 . 0-E*E ) / ( E*R) ) - { 1 . DO/E ) 

SINF=A* ( 1 . D0-E*E ) *RRD/ ( R*E*H ) 

F=DATAN2 { SINF , COSF ) 

C ECCENTRIC ANOMALY (E) 

COSEA= ( COSF+E ) / ( 1 . DO+E*COSF ) 

SINEA= ( DSQRT ( 1 . D0-E*E ) ) *SINF/ ( 1 . DO+E*COSF ) 



c 

c 


c 

c 


c 


EA=DATAN2 ( SINEA, COSEA ) 

MEAN ANOMALY (M) 

XM=EA-E* SINEA 
GEOCENTRIC ANGLE (U) 

COSU= ( X*DCOS ( OME ) +Y*DSIN ( OME ) ) /R 

SINU=( ( Y*DCOS(OME ) -X*DSIN(OME ) ) *DCOS(XI )+Z*DSIN(XI ) ) /R 
U=DATAN2 (SINU.COSU' ' ' ' 


ARGUEMENT OF PERIGEE ( omega ) 

OM=U-F 

PARTIALS OF A {semi-major axis) 

KCI (1,1)=2.0*{A**2) *X/(R**3 ) 
KCI (1,2)=2.0*(A**2) *Y/(R**3 ) 
KCI (1,3)=2.0*(A**2)*Z/(R**3) 
KCI ( 1,4 )=2. 0* (A* *2 ) *XD/GM 
KCI { 1 , 5 ) =2 . 0* ( A**2 ) * YD/GM 
KCI (1,6)=2.0*(A**2) *ZD/GM 
PARTIALS OF E (eccentricity) 

KCI ( 2 , 1 ) = { 1 . 0-E*E ) *KCI { 1 , 1 ) / ( 2 
* ( E*H*H ) 

KCI (2,2)=(1.0-E*E) *KCI ( 1 , 2 ) / { 2 
* { E*H*H ) 

KCI ( 2 , 3 ) = { 1 . 0-E*E ) *KCI ( 1 , 3 ) / ( 2 
* (E*H*H ) 

KCI (2 / 4)={1.0-E*E) *KCI ( 1 , 4 ) / { 2 
* ( E*H*H ) 


0*A*E ) - ( 1 . 0-E*E ) * { HZ*YD-HY*ZD ) / 
0*A*E) - ( 1. 0-E*E) * { HX*ZD-HZ*XD ) / 
0*A*E ) - ( 1 . 0-E*E ) * ( HY*XD-HX*YD ) / 
0*A*E)-(1.0-E*E)*( HY*Z-HZ*Y ) / 


( E*H*H ) 5 ) ( 1 * °~ E * E ) * KCI ( 1 , 5 ) / ( 2 . 0 *A*E ) - ( 1 . 0-E*E ) * ( HZ*X-HX*Z ) / 


*(e*h*h) 6, ~ (1 ' 0 ' E * E) * KCI(1 ' 6)/{2 ‘ 0 * a * e, ~ (1 * 0_e * E) * (hx * y ~ hy * x )/ 

C PARTIALS OF I (inclination) 

i I ' l \ : ! u5Iu5 *^ D ; ( hx * hx+hy *hy ) *YD ) / ( H*H*DSQRT ( HX*HX+HY*HY ) ) 

KCT ? ' R - H HY*u<7*i HX*HX+HY*HY ) * Y ) / ( H*H*DSQRT ( HX*HX+HY*HY ) ) 

O ' l _ "55*55 Z ~ ( HX * HX+ HY*HY ) *X ) / ( H*H*DSQRT ( HX*HX+HY*HY ) ) 

C PARTIMq Z * X+ 5^* HZ * Y,/(H * H * DSQRT(HX * HX+HY * HY ) ) 

c FAKllALS or OMEGA (ascending node) 

KCI ( 5 , 1 ) - ( -HX*ZD ) / ( HX*HX+HY*HY ) 

KCI ( 5 , 2 ) = ( -HY*ZD ) / ( HX*HX+HY*HY ) 

KCI ( 5 , 3 ) = ( HX*XD+HY*YD ) / ( HX*HX+HY*HY ) 

KCI ( 5 , 4 ) = ( HX*Z ) / ( HX*HX+HY*HY ) 

KCI ( 5 , 5 ) = ( HY*Z ) / ( HX*HX+HY*HY ) 

_ r,»r, mT KCI ( 5 ' 6 ) = { -HX*X-HY* Y ) / ( HX*HX+HY*HY ) 

C PARTIALS of M (mean anomaly) 

*SINF ) t w ! * KCI ( 1 • 1 )/( R*E ) -2 . *A* (C0SF*RRD/H- 

* WH+ »??? n ' IIC w o ( 1 ‘ 0_E E ) / ( R*E*H ) ) * ( C0SF*RRD/H ) * ( HZ*YD-HY*ZD 
) /H+ ( A* ( 1 . 0-E*E ) / ( R*E*H ) ) *C0SF*XD+ ( SINF*X/ ( R*R*E ) ) 

*SXNF!iKci?? E 2t/S 0 fn^ tD n H ;?Jw > * KCI(1 ' 2)/(R * E) - 2 -* A * (COSF '* RRD /H- 

* ; KCI ^ 2 ' ^ )/ R ~ ( A* ( 1 • 0-E*E )/(R*E*H ) ) * (C0SF*RRD/H ) * (HX*ZD-HZ*XD 
)/H+(A*(1.0-E*E)/(R*E*H) ) *C0SF*YD+ ( SINF*Y/ ( R*R*E ) ) 

♦SINF ) /R°f n^?°n H p S ^ w } * KCI ( 1 ' 3 ) / ( R * E ) -2 . *A* ( C0SF*RRD/H- 
*f /5+ a n ' HC w A ( 1 • 0_E E ) / ( R*E*H ) ) * ( C0SF*RRD/H ) * { HY*XD-HX*YD 

)/H+( A* ( 1 . 0-E*E)/(R*E*H) ) *C0SF*ZD+ ( SINF*Z/ (R*R*E ) ) 

*SINF S *KC-H 2 E i t w * * KCI (1 ' 4)/(R * E >“ 2 * * A * ( COSF*RRD/H- 

*/H + (|. ( l! 0 -E*e(/ ( R* E iJ i ?4f&T E * H ’ ’ * (COSF.RRD/H ) * ( HY*Z-HZ*Y ) 

*smf ! 1 ! * KC1 < ■ 1 > 1 5 >/< r * e > -2 • *a* < cosf.rwh- 

^ ( 2 / 5 ) /R- ( A* ( 1 . 0-E^E ) / ( R*E*H ) ) * ( C0SF*RRD/H ) * { HZ*X-HX*Z ) 
*/H+ ( A* { 1 . 0-E*E ) / ( R*E*H ) ) *COSF*Y 7 ; 

F6= ( 1 . 0-E*E ) * ( COSF*RRD/H-SINF ) *KCI (1,6 )/( R*E ) -2 . *A* ( COSF*RRD/H- 



•SINF ) *KCI ( 2 , 6 )/R- < A* ( 1 . 0-E*B ) / ( R*E*H )) * (COSF*RED/H ) * ( HX*Y-HY*X) 
*/H+(A*(l. 0-E*E ) / ( R*E*H > > *C0SF*2 


E4=( (DSQRT( 1 . t ' n+F*COSF) ) * (F5-ilNF*KCI ( 2, 5 ) /( 1 • 0 - E*E) , 

KCI (6 1 ) = ( 1 . 0-E*COSEA ) *E1 -SINEA*KCI (2,1) 

KCI ( 6 # 2 ) = ( 1 • 0-E*COSEA) *E2-SINEA*KCI (2,2) 

KCI 6 3 - 1 . 0-E*COSEA ) *E3-SINEA*KCI (2,3 
KCI 6 i 4 = ( 1 . 0-E*COSEA ) *E4-SINEA*KCI (2,4 
KCI (6 5 ) =( 1 . 0-E*COSEA) *E5-SINEA*KCI (2,5) 

KCI ( 6 ! 6 ) = ( 1 • 0-E* COSEA ) *E6-SINEA*KCI (2,6) 

PARTIALS of omega . * COS U-X*DCOS (OME) *DCOS (XI ) * 

siSSfeass;., 

:f + {-^™(0°M^ : Dcjs,X i!*^S^^IHiicOS ( OME,.DOO S( XI,. 

“ tSS(OMB) *DC0S(X1 [ ) ^““{SflSSSiiDcoa.ao ) *DCOS (XI . . 

iiil isSSEii. 


*) ) 

KCI (4,1) =U1-F1 
KCI (4,2) =U2-F2 
KCI (4,3) =U3-F3 
KCI (4,4) =U4-F4 
KCI (4,5) =U5-F5 
KCI (4,6) =U6-F6 
RETURN 
END 



***INPUT DATA: AO , EO , XIO ,OMO , OMEO , XMO , FO ,EAO ,XO , YO , ZO , XDO , YDO ^^**** 


* 17570925.77037010 
*-0.881968672701392 

* 1.929121467522410 
*-7615946.972994897 

* 2396.747150406745 


0.568848843824957 

3.124086383674180 

1.294798081992500 

-10889534.81787064 

-4319.688274669979 


0.541052068118242* 

0.747478091842417* 

* 

6622196.773840833* 

2569.923171671754* 

********************** 


************* 


**** ELEMENTS OF THE SATELLITE STATE VECTOR (XS,YS,ZS,XSD,YSD,ZSD) ***** 
-7615946.97299500000 -10889534.81787060000 6622196.77384083000 

2396.74715040671 -4319.68827467000 2569.92317167176 


LINEAR KEPLER TRANSFORMATION MATRIX KC ******** 


******** ELEMENTS OF THE 
-0.4334402792730 
12744314.0219510000000 
-0.6197473576625 
-6587846.0049455100000 
0.3768837715431 
3823722.2809170100000 
-0.0000682020737 
5026 . 1010036825600 
0.0001229214764 
2031.2430250534000 
-0.0000731299877 
-1273.1729165956700 
********************** 


28607803.1683634000000 

10889534.8178706000000 

-20941688.6224014000000 

-7615946.9729950000000 

12280202.6842647000000 

0.0000000000000 

4895.4299974719600 

4319.6882746700000 

3878.7401494792500 

2396.7471504067100 

-2381.7166470278200 

0.0000000000000 


115924.0427484690000 

8841893.9377193300000 

6621182.0491344800000 

-15935859.3843162000000 

11021186.2187911000000 

9480761.5013344300000 

44.9874707423364 

3421.8010333312000 

2569.5293802118000 

4892.6051644543400 

4277.0704058304800 

-2975.3148024795100 




****** KEFLERIAN ORBITAL ELEMENTS (ASfES^ISjOMS^MESfXMSfFS^AS) ****** 
1.757092577037010D+007 0.568848843824957 0.541052068118242 

-0.881968672701392 3.124086383674180 0.747478091842417 

1.929121467522410 1.294798081992500 


************ INVERSE LIN] 
-1.43685580235904000 
3712.81862519266000000 
0.00000003821214436 
-0.00009454042695639 
0.00000000032380779 
0.00000097249246146 
-0.00000008348546989 
0.00049142428946580 
-0.00000000126905576 
0.00000327011214863 
0.00000004091773235 
-0.00022446651641490 


R KEPLER TRANSFORMATION 
-2.05446431593198000 
6691.66085312722000000 
-0.00000003442479720 
-0.00011261293567051 
0.00000001849478572 
0.00005554541988079 
-0.00000000487797231 
-0.00052188742069056 
-0.00000007248409427 
0.00018677754281016 
0.00000010672721602 
0.00055778897052638 


MATRIX KCI *********** 
1.24937081266399000 
3981.08687246287000000 
0.00000002027940900 
0.00006864871363555 
0.00000003078520965 
0.00009245726995033 
0.00000014456563439 
-0.00005433900188108 
-0.00000012065227854 
0.00031089767137094 
-0.00000006454873928 
-0.00033274106616196 


FULFILMENT OF KC*KCI=I 

*********************************************************************** 

1.00000000000000000000 0.00000000000000591861 - 0.00000000000000149414 

- 0.00000000000725564053 - 0.00000000000042099657 0.00000000000037148062 

*************************************** 

- 0.00000000000001013875 1.00000000000000000000 0.00000000000000232355 

0.00000000000916000609 - 0.00000000001096189806 - 0.00000000000075761619 

*************************************** 

0.00000000000000400881 0.00000000000000253020 1.00000000000000000000 

- 0.00000000000583844084 - 0.00000000000010880186 - 0.00000000001283706474 

*************************************** 

0.00000000000000000072 - 0.00000000000000000120 0.00000000000000000070 

1.00000000000000000000 0.00000000000000301040 0.00000000000000069714 

*************************************** 

- 0.00000000000000000062 0.00000000000000000096 0.00000000000000000043 

- 0.00000000000000469340 1.00000000000000000000 - 0.00000000000000116671 

*************************************** 

0.00000000000000000041 0.00000000000000000052 0.00000000000000000138 

0.00000000000000055576 - 0.00000000000000099508 1.00000000000000000000 

*********************************************************************** 


************* 


FULFILMENT OF KCI*KC=I 

*********************************************************************** 

1.00000000000000000000 0.00000137576353154145 0.00000000985528458841 

0.00000062681783674634 0.00000055624150263611 0.00000072245620685862 

*************************************** 

0.00000000000000000000 1.00000000000001000000 0.00000000000000018922 

0.00000000000000450163 0.00000000000000328440 - 0.00000000000000115492 

*************************************** 

0.00000000000000000000 - 0.00000000000000293900 1.00000000000000000000 

- 0.00000000000000185891 - 0.00000000000000358525 - 0.00000000000000122450 

*************************************** 

0.00000000000000000000 0.00000000000007482006 - 0.00000000000000375618 

1.00000000000004000000 0.00000000000003101946 0.00000000000004706265 

*************************************** 

0.00000000000000000000 0.00000000000000450378 0.00000000000000818540 

0.00000000000000009552 1.00000000000000000000 0.00000000000000027067 

********. V ****************************** 

0.00000000000000000000 - 0.00000000000008301422 - 0.00000000000000178362 
- 0.00000000000003955419 - 0.00000000000003490144 0.99999999999996000000 

*********************************************************************** 


